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Abstract 


This  work  develops  an  adaptive  concurrent  multi-level  computational  model  for  multi-scale  analysis  of 
composite  structures  undergoing  damage  initiation  and  growth  due  to  microstructural  damage  induced  by 
debonding  at  the  fiber-matrix  interface.  The  model  combines  macroscopic  computations  using  a  continuum 
damage  model  developed  in  a  preceding  paper  [75]  with  explicit  micromechanical  computations  of  stresses 
and  strain,  including  explicit  debonding  at  the  fiber-matrix  interface.  The  macroscopic  computations  are 
done  by  conventional  FEM  models  while  the  Voronoi  cell  FEM  is  used  for  micromechanical  analysis.  Three 
hierarchical  levels  of  different  resolution  adaptively  evolve  in  this  to  improve  the  accuracy  of  solutions  by 
reducing  modeling  and  discretization  errors.  They  levels  include:  (a)  level-0  of  pure  macroscopic  analy¬ 
sis  using  a  continuum  damage  mechanics  (CDM)  model;  (b)  level-1  of  asymptotic  homogenization  based 
macroscopic-microscopic  RVE  modeling  to  monitor  the  breakdown  of  continuum  laws  and  signal  the  need 
for  microscopic  analyses;  and  (c)  level-2  regions  of  pure  micromechanical  modeling  with  explicit  depiction  of 
the  local  microstructure.  Two  numerical  examples  are  solved  to  demonstrate  the  effectiveness  and  accuracy 
of  the  multi-scale  model.  A  double  lap  bonded  composite  joint  is  modeled  for  demonstrating  the  model’s 
capability  in  handling  large  structural  problems. 

For  micromechanical  analysis,  an  extended  Voronoi  cell  finite  element  model  (X-VCFEM)  is  developed 
for  modeling  multiple  cohesive  crack  propagation  in  brittle  materials.  The  cracks  are  modeled  by  a  cohesive 
zone  model  and  their  incremental  directions  and  growth  lengths  are  determined  in  terms  of  the  cohesive 
energy  near  the  crack  tip.  Extension  to  VCFEM  is  achieved  through  enhancements  in  stress  functions  in 
the  assumed  stress  hybrid  formulation.  In  addition  to  polynomial  terms,  the  stress  functions  include  branch 
functions  in  conjunction  with  level  set  methods,  and  multi-resolution  wavelet  functions  in  the  vicinity  of 
crack  tips.  Comparison  of  X-VCFEM  simulation  results  with  results  in  literature  for  several  fracture  me¬ 
chanics  problems  validates  the  effectiveness  of  X-VCFEM.  Effect  of  stereographic  features  such  as  size  and 
distribution  of  heterogeneities  on  damage  evolution  in  random  microstructures  are  also  discussed.  In  order 


to  study  the  interaction  between  interface  debonding  and  cohesive  matrix  cracking,  a  criterion  based  on 
cohesive  models  is  proposed  to  assess  the  crack  penetrating  into  matrix  from  the  interface  and  is  validated 
by  numerical  examples. 

Next  the  extended  Voronoi  cell  finite  element  model(X-VCFE]V[)  has  been  developed  for  modeling  in¬ 
terfacial  debonding  with  arbitrary  matrix  cohesive  cracking  in  fiber-reinforced  composites.  To  describe  the 
onset  and  growth  of  damage  along  the  fiber-matrix  interface,  normal  and  tangential  cohesive  zone  models 
are  coupled  into  VCFEM.  It  is  shown  that  the  initiation  and  especially  propagation  of  debonding  depends 
not  only  on  the  total  cohesive  energy,  but  also  on  the  shape  of  the  traction-displacement  curve.  The  model 
is  also  used  to  study  the  influence  of  various  local  morphological  parameters  on  damage  evolution  by  inter¬ 
facial  debonding.  A  special  function  of  various  geometric  parameters  is  developed  to  predict  the  location  of 
debonding  in  microstructures  with  varying  morphology. 

Finally,  a  Voronoi  cell  finite  element  model  is  also  developed  for  transient  elastodynamic  analysis  in  time 
domain  is  developed.  In  the  present  formulation,  the  inertia  field  is  approximated  in  terms  of  stresses  so 
as  to  satisfy  the  equilibrium  equation  a-priori.  The  weak  forms  of  kinematics  and  traction  reciprocity  are 
obtained  by  minimization  of  the  complementary  variational  principle.  Stress  wave  is  a  local  disturbance  that 
propagates  through  the  material,  resulting  in  high  stress  gradients  near  the  wave  front.  Therefore,  localization 
and  multi-resolution  properties  of  the  wavelet  functions  are  exploited  to  enhance  the  computational  efficiency 
by  enriching  the  stress  functions  only  locally  near  the  wave  front.  The  enrichment  is  carried  out  adaptively 
by  employing  posteriori  local  error  estimators  that  determine  the  required  translation  and  dilation  of  the 
wavelet  functions  at  each  time  step.  The  accuracy  and  computational  efficiency  of  the  proposed  method  is 
demonstrated  through  comparison  with  analytical  solutions  and  conventional  FEM  packages. 


2 


Chapter  1 


Relevant  Information 

1.1  Personnel  Supported 

1.  Somnath  Ghosh,  PI 

2.  Jayesh  Jain,  Ph.D.  Student,  Graduate  Research  Associate,  100% 

3.  Shanhu  Li,  Ph.D.  Student,  Graduate  Research  Associate,  (partial  support) 

Now  at  ALGOR  Corporation. 

4.  Shriram  Swaminathan,  M.S.  Student,  Graduate  Research  Associate,  (partial  support) 

Now  at  ABAQUS 

1.2  Completed  Ph.D.  Dissertations  &  M.S.  Thesis 

1.  J.  Jain,  Ph.D.  candidate.  Dissertation  Topic:  Multiple  Scale  Modeling  of  Damage  in  Composite  Mate¬ 
rials  under  Dynamic  Conditions 

2.  S.  Li,  Ph.D.  2005,  Dissertation  Title:  Extended  Voronoi  Cell  Finite  Element  Model  for  Damage  in 
Brittle  Matrix  Composites 


1 


3.  S.  Swaminathan  M.S.  2005,  Thesis  Title:  Statistically  Equivalent  Representative  Volume  Elements  for 
Composite  Microstructures 

1.3  Refereed  Journal  Publications 

1.  J.  R.  Jain  and  S.  Ghosh,  A  3D  continuum  damage  mechanics  model  for  fiber  reinforced  composites 
with  interfacial  damage,  submitted  for  publication. 

2.  J.  R.  Jain  and  S.  Ghosh,  Multi-resolution  wavelet  enriched  Voronoi  cell  finite  element  model  for  elastic 
wave  propagation  in  heterogeneous  solids,  in  preparation. 

3.  S.  Li  and  S.  Ghosh,  Modeling  interfacial  debonding  and  matrix  cracking  in  fiber  reinforced  composites 
by  the  extended  Voronoi  cell  EEM,  submitted  for  publication. 

4.  S.  Ghosh,  J.  Bai  and  P.  Raghavan,  Concurrent  multi-level  model  for  damage  evolution  in  microstruc- 
turally  debonding  composites.  Mechanics  of  Materials,  (in  press). 

5.  S.  Li  and  S.  Ghosh,  Multiple  cohesive  crack  growth  in  brittle  materials  by  the  extended  Voronoi 
cell  finite  element  model.  International  Journal  of  Eracture,  (in  press).  Published  online  2006:  DOI: 
10:1007/sl0704-006-9000-2. 

6.  S.  Li  and  S.  Ghosh,  Extended  Voronoi  cell  finite  element  model  for  multiple  cohesive  crack  propagation 
in  brittle  materials.  International  Journal  for  Numerical  Methods  in  Engineering  Vol.  65,  pp.  1028- 
1067,  2006. 

7.  P.  Raghavan  and  S.  Ghosh,  A  continuum  damage  mechanics  model  for  unidirectional  composites  un¬ 
dergoing  interfacial  debonding.  Mechanics  of  Materials,  Vol.  37,  No.  9,  pp.  955-979,  2005. 

8.  S.  Swaminathan,  S.  Ghosh  and  N.J.  Pagano,  Statistically  equivalent  representative  volume  elements 
for  composite  microstructures.  Part  I:  Without  damage.  Journal  of  Composite  Materials  Vol.  40,  No. 
7,  pp.  583-604,  2006. 


2 


9.  S.  Swaminathan,  S.  Ghosh  and  N.J.  Pagano,  Statistically  equivalent  representative  volume  elements 


for  composite  microstructures,  Part  II:  With  evolving  damage,  Journal  of  Composite  Materials,  Vol. 
40,  No.  7,  pp.  605-621,  2006. 

10.  S.  Ghosh  and  P.  Raghavan,  Multi-scale  model  for  damage  analysis  in  fiber  reinforced  composites  with 
interfacial  debonding  materials.  International.  Journal  of  Multiscale  Computational  Engineering,  Vol. 2, 
No.  4,  pp.  621-645,  2004. 

11.  S.  Ghosh  and  S.  Moorthy,  Three  dimensional  Voronoi  cell  finite  element  model  for  microstructures 
with  ellipsoidal  heterogeneities”,  Vol.  34,  pp.  510-531,  2004. 


1.4  Proceeding  Publications 

1.  S.  Ghosh  and  S.  Li,  A  wavelet  based  Voronoi  cell  FEM  for  crack  evolution  in  composite  materials. 
Proceedings  of  the  11th  International  Congress  on  Fracture,  A.  Carpinteri  (eds.),  Turin,  Italy,  pp. 
4905-5001,  2005. 

2.  S.  Ghosh,  Multiple  scale  modeling  for  deformation  and  failure  of  heterogeneous  materials.  Proceedings 
of  the  International  Congress  on  Computational  Mechanics  and  Simulation,  NCR  Iyengar  and  A. 
Kumar  (eds.),  IIT  Kanpur  Press,  pp.  43-56,  2004. 

3.  S.  Ghosh  Modeling  at  the  interface  of  mechanics  and  materials  for  composite  and  polycrystalline 
materials  Proceedings  ICTACEM  2004;  Third  International  Congress  on  Theoretical,  Applied,  Com¬ 
putational  and  Experimental  Mechanics,  S.  Bhattacharyya  and  S.  Ghosh  (eds.),  IIT  Kharagpur  Press, 
pp.  10-11,  2004. 

4.  P.  Raghavan,  J.  Bai  and  S.  Ghosh,  Multi-scale  model  for  damage  analysis  in  fiber-reinforced  composites 
with  debonding.  Materials  Processing  and  Design:  Modeling,  Simulation  and  Application  Proceedings 
of  NUMIFORM,  S.  Ghosh,  J.  Castro  and  J.K.  Lee  (eds.),  AIP  Publishers,  pp.  1911-1917,  2004. 


3 


1.5 


Keynote  Lectures 


1.  S.  Ghosh,  Computational  multi-scale  models  for  structure-material  interaction,  AFOSR  Workshop  on 
Multiscale  Modeling  of  Carbon  Fiber  Reinforced  Composites,  Long  Beach,  CA,  May,  2006. 

2.  S.  Ghosh,  Computational  multi-scale  models  for  structure-material  interaction.  Mechanics  of  Materials 
Workshop,  Mathematisches  Forschungsinstitut,  Oberwolfach,  Germany,  January,  2006. 

3.  S.  Ghosh,  Concurrent  multi-scale  computational  models  for  structure-material  interaction,  Ohio  Su¬ 
percomputer  Center  on  Materials  for  National  Security,  Santa  FE,  New  Mexico,  May,  2005 

4.  S.  Ghosh  Multiple  Scale  Modeling  for  Deformation  and  Failure  of  Heterogeneous  Materials,  Interna¬ 
tional  Congress  on  Computational  Mechanics  and  Simulation,  Indian  Institute  of  Technology,  Kanpur, 
December  2004. 

5.  S.  Ghosh,  P.  Raghavan  and  S.  Li,  Multi-Level  Computational  Models  For  Multiple  Scale  Analysis  of 
Composite  Materials  International  Workshops  on  Advances  in  Computational  Mechanics:  IWACOM, 
Hosei  University,  Tokyo,  Japan,  November  2004. 

6.  S.  Ghosh  and  P.  Raghavan,  ‘  Multi-Level  Models  For  Damage  Analysis  In  Composite  Materials,  7th 
US  National  Congress  of  Computational  Mechanics,  Albuquerque,  NM,  July  2003. 

7.  S.  Ghosh  and  P.  Raghavan,  ‘Adaptive  Multi-Level  Models  For  Damage  Analysis  In  Composite  Mate¬ 
rials,  International  Conference  on  Computational  Engineering  Science  ICES  2003,  Corfu,  Greece,  July 
2003. 

8.  S.  Ghosh,  ‘  Multilevel  Models  for  Multiple  Scale  Analysis  of  Composite  Materials’,  Annual  Symposium 
in  the  Computational  Science  and  Engineering  Program,  The  University  of  Illinois,  Urbana-Champaign, 
April  2003 


4 


1.6  Invited  Lectures  at  Conferences 

1.  S.  Ghosh,  J.  Bai  and  P.  Raghavan,  Concurrent  multi-level  model  for  damage  evolution  in  microstruc- 
turally  debonding  composites,  7th  World  Congress  of  Computational  Mechanics,  Los  Angeles,  CA, 
July,  2006. 

2.  S.  Li  and  S.  Ghosh,  Extended  Voronoi  cell  finite  element  model  for  damage  evolution  in  composite 
materials,  7th  World  Congress  of  Computational  Mechanics,  Los  Angeles,  CA,  July,  2006. 

3.  S.  Li  and  S.  Ghosh,  Extended  Voronoi  cell  finite  element  model  (X-VCEEM)  for  multiple  cohesive  crack 
propagation  in  brittle  materials  ,8th  US  National  Congress  on  Computational  Mechanics,  Austin,  TX, 
July,  2005. 

4.  S.  Ghosh  and  S.  Li,  An  Extended  Voronoi  cell  EEM  (X-VCEEM)  for  multiple  cohesive  crack  evolution 
in  composite  materials,  11th  International  Congress  on  Eracture,  Turin,  Italy,  March,  2005. 

5.  S.  Ghosh,  P.  Raghavan,  J.  Bai  and  S.  Li,  Multi-Level  computational  models  for  multiple  scale  analysis 
of  composite  materials  ,  International  Mechanical  Engineering  Conference  and  Exposition,  IMECE 
2004,  Anaheim,  California,  November  2004. 

6.  S.  Ghosh,  Multi-Scale  model  for  damage  analysis  in  fiber  reinforced  composites  with  debonding  Mul¬ 
tiscale  Modeling  Conference:  MMM-2,  University  of  California,  Los  Angeles,  October  2004. 

7.  S.  Li  and  S.  Ghosh,  A  wavelet  based  Voronoi  Cell  EEM  for  crack  evolution  in  composite  materials  6th 
World  Congress  of  Computational  Mechanics:  WCCM,  Beijing,  China,  September  2004. 

8.  P.  Raghavan,  J.  Bai  and  S.  Ghosh,  Multi-scale  model  for  damage  analysis  in  fiber-reinforced  composites 
with  debonding,  NUMIEORM  2004:  8th  International  Conference  on  Numerical  Methods  in  Industrial 
Eorming  Processes,  The  Ohio  State  University,  Columbus  Ohio,  June  2004. 

9.  S.  Ghosh  and  P.  Raghavan,  ‘  Multilevel  models  for  multiscale  damage  analysis  in  composite  materials 
’,  2nd  MIT  Conference,  Cambridge  MA,  June  2003. 
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1.7 


Honors  and  Awards 


1.  S.  Ghosh:  Fellow,  ASM  International,  The  Materials  Information  Society,  2006 

2.  S.  Li:  (Ph.D.  student  and  GRA):  Robert  J.  Melosh  Medal  Competition  Finalist,  2006 

3.  S.  Manchiraju:  (Ph.D.  student  and  GRA):  Winner  of  Materials  Modeling  Student  Competition,  at  the 
7th  World  Congress  of  Computational  Mechanics  at  Los  Angeles,  in  2006 

4.  S.  Ghosh:  John  B.  Nordholt  Professorship,  College  of  Engineering,  The  Ohio  State  University,  2004 

5.  S.  Ghosh:  Lumley  Interdisciplinary  Research  Award,  College  of  Engineering,  The  Ohio  State  University, 
2004 

6.  S.  Ghosh:  Chairman  of  NUMIEORM  2004:  Numerical  Methods  in  Industrial  Eorming  Processes’,  June 
2004. 

7.  S.  Ghosh:  Lumley  Eaculty  Research  Award,  College  of  Engineering,  The  Ohio  State  University,  2003 

8.  P.  Raghavan:  Robert  J.  Melosh  Medal  Winner,  for  best  student  paper  in  Einite  Element  Analysis, 
2003 

9.  Member  of  Executive  Council,  U.S.  Association  of  Computational  Mechanics,  2002- 

1.8  Collaborations  with  Army  Research  Laboratory 

We  have  collaborated  with  Dr.  Peter  Chung  and  Dr.  Raju  Namburu,  at  the  Army  Research  Laboratory. 
The  PI  has  visited  ARL  several  times  to  present  seminar  and  discuss  progress  with  researchers  in  the 
Computational  and  Informational  Sciences  Directorate  at  ARL.  Dr.  Namburu  and  Chung  have  also  visited 
OSU  twice  to  monitor  the  breadth  and  depth  of  progress  at  the  Pis  laboratory. 
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1.9  Technology  Transfer  to  Army  Research  Laboratory 


Technology  transfer  of  some  of  the  codes  developed  in  the  Pis  research  laboratory  to  ARL  platforms  has 
taken  place  in  consultation  with  Dr.  Peter  Chung  and  Dr.  Raju  Namburu  in  the  Computational  and  Infor¬ 
mational  Sciences  Directorate  at  ARL.  The  PI  visited  has  ARL  several  times  in  the  last  few  years  and  Dr. 
Namburu  and  Dr.  Chung  have  also  visited  the  Computational  Mechanics  Laboratory  at  OSU  several  times. 
The  recent  visit  by  the  PI  to  ARL  in  2006  was  on  July  14,  2006.  Significant  exchange  of  ideas  and  future 
plans  took  place  during  this  visit. 

The  PI  has  transferred  major  codes  to  the  ARL  platforms  for  use  by  ARL  researchers.  A  special  person 
with  security  clearance  and  authorized  to  work  on  US  government  systems  has  been  recruited  from  the  Ohio 
Supercomputer  Center  to  facilitate  the  technology  transfer  process.  The  first  code  that  has  been  transferred 
is  the  crystal  plasticity  model.  Initial  efforts  were  identified  to  provide  access  to  the  crystal  plasticity  model 
through  a  user  defined  material  model  in  LS-DYNA.  This  required  the  proper  LSTC  software,  a  special  ver¬ 
sion  of  LS-DYNA  that  allows  user  subroutines  to  be  compiled  in,  to  be  verified  on  the  target  ARL  systems. 
These  systems  are  identified  to  be  the  SGI  Origin,  the  IBM  P3  system  and  the  Linux  cluster.  Only  the 
Origin  had  the  software  installed,  so  requests  were  entered  through  the  ARL  support  staff  and  the  software 
is  now  installed  on  these  systems.  Next,  the  user  subroutines  were  moved  to  the  three  ARL  systems  and 
compilation  and  numerical  accuracy  needed  to  be  verified.  Compilation  scripts  have  been  customized  for 
each  system  and  a  validation  model  was  run  so  that  numerical  results  could  be  compared.  Compilation 
was  verified  on  all  three  machines  and  numerical  results  have  been  verified  on  the  Origin  and  IBM  systems. 
Documentation  is  also  provided  for  users  to  be  able  to  reference  a  README  file  on  each. 

The  second  code  and  model  that  has  been  ported  to  ARL  systems  is  an  earlier  version  of  the  CDM 
model.  The  Eortran  90  code  has  been  ported  to  JVN  Linux  cluster  at  ARL.  We  have  created  Make  files 
to  build  the  model  and  have  created  batch  scripts  documenting  example  usage.  The  User  Manual  for  this 
code  documents  important  control  variables,  input  file  structure  and  provides  an  example  case  to  document 
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analysis  workflow. 
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Chapter  2 


Introduction 


Analysis  of  composite  materials  with  microstructural  heterogeneities  is  conventionally  done  with  macroscopic 
properties  obtained  by  homogenizing  response  functions  in  the  representative  volume  element  (RVE)  from 
microscopic  analyses  at  smaller  length  scales.  While  these  “bottom-up”  homogenization  models  are  efficient 
and  can  reasonably  predict  macroscopic  or  averaged  behavior,  such  as  stiffness  or  strength,  they  have  limited 
predictive  capabilities  with  problems  involving  localization,  failure  or  instability.  Assumptions  of  macroscopic 
uniformity  and  RVE  periodicity,  the  two  basic  requirements  of  homogenization,  break  down  under  these 
circumstances.  The  uniformity  assumption  ceases  to  hold  in  critical  regions  of  high  local  solution  gradients, 
such  as  near  free  edges,  interfaces,  material  discontinuities  or  evolving  damage.  RVE  periodicity,  on  the  other 
hand,  is  unrealistic  for  non-uniform  microstructures,  e.g.  in  the  presence  of  clustering  of  heterogeneities  or 
microscopic  damage.  Even  with  a  uniform  phase  distribution  in  the  microstructure,  the  evolution  of  localized 
stresses,  strains  or  damage  path  can  violate  the  periodicity  conditions.  Problems  like  this  have  been  effectively 
tackled  by  multi-scale  modeling  methods  e.g.  in  [72,  29,  44,  67,  66,  82,  81,  74,  73,  94,  108,  92].  Multi-scale 
analyses  methods  can  be  broadly  classified  into  two  classes.  The  first  is  known  as  ’’hierarchical  models” 
[29,  44,  94,  92]  in  which  information  is  passed  from  lower  to  higher  scales,  usually  in  the  form  of  material 
properties.  The  hierarchical  homogenization  models  assume  periodic  representative  volume  elements  (RVE) 
in  the  microstructure  and  uniformity  of  macroscopic  field  variables.  The  second  class,  known  as  “concurrent 
methods”  [81,  67,  66,  82,  74,  73,  108],  implement  sub-structuring  and  simultaneously  solve  different  models 
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at  regions  with  different  resolutions  or  scales. 

Two-way  coupling  of  scales  enabled  in  the  concurrent  methods  is  suitable  for  problems  involving  local¬ 
ization,  damage  and  failure.  Macroscopic  analysis,  using  bottom-up  homogenization  in  regions  of  relatively 
benign  deformation,  enhances  the  efficiency  of  the  computational  analysis.  As  a  matter  of  fact,  it  would  be 
impossible  to  analyze  large  structural  regions  without  the  advantage  of  a  continuum  model  based  macro¬ 
scopic  analysis.  On  the  other  hand,  the  top-down  localization  process  cascading  down  to  the  microstructure 
in  critical  regions  of  localized  damage  or  instability  for  pure  microscopic  analysis,  is  necessary  for  accu¬ 
rately  predicting  the  damage  path.  These  microscopic  computations,  depicting  the  real  microstructure  are 
often  complex  and  computationally  prohibitive.  Hence,  a  concurrent  setting  makes  such  analyses  feasible, 
provided  the  ” zoom-in”  regions  are  kept  to  a  minimum.  The  adaptive  multi-level  models,  promoted  in 
[67,  66,  82,  74,  73,  108],  are  attempts  to  achieve  this  objective,  with  the  adaptivity  motivated  from  physi¬ 
cal  and  mathematical  perspectives.  However,  there  is  a  paucity  of  such  studies  in  the  literature  involving 
material  nonlinearity  and  evolving  microstructural  damage.  In  their  previous  studies,  Ghosh  and  coworkers 
have  proposed  adaptive  multi-level  analysis  using  the  microstructural  Voronoi  cell  FEM  model  for  modeling 
elastic-plastic  composites  with  particle  cracking  and  porosities  in  [81],  and  for  elastic  composites  with  free 
edges  and  stress  singularities  in  [74,  73]. 

In  this  work,  we  have  derived  and  computationally  modeled  an  anisotropic  continuum  damage  mechanics 
(CDM)  model  for  unidirectional  fiber-reinforced  composites  undergoing  interfacial  debonding  from  by  using 
homogenization  theory.  The  CDM  model  homogenizes  the  damage  incurred  through  initiation  and  growth 
of  interfacial  debonding  in  a  microstructural  RVE  with  nonuniform  distribution  of  fibers.  Additionally, 
arbitrary  loading  conditions  are  also  effectively  handled  by  this  model.  The  CDM  model  is  then  used  in  an 
adaptive  concurrent  multi-level  computational  model  to  analyze  multi-scale  evolution  of  damage.  Damage 
by  fiber-matrix  interface  debonding,  is  explicitly  modeled  over  extended  microstructural  regions  at  critical 
locations  [35,  53].  The  adaptive  model  addresses  issues  of  efficiency  and  accuracy  through  considerations  of 
physically-based  modeling  errors. 

The  adaptive  multi-level  model  consists  of  three  levels  of  hierarchy  viz.  level-0,  level- 1  and  level-2),  which 
evolve  in  sequence.  The  continuum  damage  model  developed  in  [75]  is  used  for  level-0  computations.  The 
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level-1  domain  is  used  as  a  ‘swing  region’  to  establish  criteria  for  switching  from  macroscopic  to  microscopic 
calculations.  Physical  criteria  involving  variables  at  the  macroscopic  and  microstructural  RVE  levels,  trigger 
switching  from  pure  macroscopic  to  pure  microscopic  calculations,  i.e.  the  level  —  0  — >•  level  —  1  — >•  level  —  2. 
A  transition  layer  is  placed  between  the  level  —  1  and  microscopic  level  —  2  domains  for  smooth  transition 
from  one  scale  to  the  next. 

An  important  damage  phenomenon  in  composite  microstructures  is  crack  propagation  in  brittle  matrix. 
Numerical  analysis  and  simulation  of  the  growth  of  multiple  cracks  in  materials  is  a  challenging  enterprise 
due  to  morphological  and  constitutive  complexities  that  govern  its  growth.  Even  a  very  high  density  mesh 
cannot  overcome  pathological  mesh  dependence  near  the  crack  tips  and  avoid  biasing  the  direction  of  crack 
propagation.  The  difficulties  aggravate  in  the  presence  of  multiple  cracks,  due  to  their  interaction  with  each 
other.  Various  methods  have  been  proposed  for  improving  the  effectiveness  of  computational  methods  in 
modeling  cracks.  While  most  of  these  analyses  are  limited  to  stationary  cracks,  it  is  only  recently  that  effec¬ 
tive  methods  of  analysis  of  crack  propagation  are  being  proposed.  With  increasing  power  of  computational 
modeling  and  hardware,  the  cohesive  zone  models  [63,  64,  65,  97,  30,  33,  39,  68]  have  emerged  as  important 
tools  for  modeling  crack  propagation  in  homogeneous  and  heterogeneous  materials.  In  these  models,  inter¬ 
faces  of  similar  and  dissimilar  materials  are  treated  as  zero  thickness  non-linear  springs.  Interfacial  traction 
is  specified  as  nonlinear  functions  of  tangential  and  normal  separations  across  the  interface  to  manifest  crack 
evolution.  These  models  have  been  used  to  simulate  crack  growth  between  elements  in  [13,  103,  39],  by  lacing 
the  interface  between  contiguous  elements  with  cohesive  springs.  The  use  of  a  highly  refined  computational 
mesh,  especially  near  the  crack  tip  is  still  a  requirement,  even  though  the  effect  is  mitigated  due  to  the 
finite  crack  tip  stress  with  this  model.  Alternatively,  intra-element  enrichment  approaches,  based  on  the 
incorporation  of  embedded  discontinuities  in  displacement  or  strain  fields  have  been  proposed  ([48]),  which 
eliminates  mesh  dependent  prediction  of  the  evolving  crack  path,  and  hence  the  need  for  remeshing.  The 
extended  EEM  or  X-EEM  [8,  7,  9,  10,  24,  59,  60]  is  a  powerful  recent  addition  to  this  family  of  intra-element 
enrichment.  Cohesive  crack  propagation  has  been  modeled  in  this  work  by  using  the  partition  of  unity 
concept  to  incorporate  local  enrichment  functions  that  allows  the  preservation  of  the  general  displacement 
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based  FEM  formalism. 


Stress-based  finite  element  methods  have  had  considerable  success  when  stress  fields  are  of  interest  in  the 
analysis  [96,  95].  Within  this  general  formalism,  the  Voronoi  cell  finite  element  method  (VCFEM)  has  been 
developed  in  [35,  61,  80,  36,  79,  78,  53]  for  micromechanical  analysis  of  arbitrary  heterogeneous  microstruc¬ 
tures.  The  method  can  effectively  overcome  requirements  of  large  degrees  of  freedom  in  conventional  finite 
element  models.  Morphological  arbitrariness  in  dispersions,  shapes  and  sizes  of  heterogeneities,  as  seen  in 
real  micrographs  are  readily  modeled  by  this  method.  The  VCEE  model  naturally  evolves  by  tessellation 
of  the  microstructure  into  a  network  of  multi-sided  Voronoi  polygons.  Each  Voronoi  cell  with  embedded 
heterogeneities  (particle,  fiber,  void,  crack  etc.)  represents  the  region  of  contiguity  for  the  heterogeneity, 
and  is  treated  as  an  element  in  VCEEM.  VCEEM  elements  are  considerably  larger  than  conventional  EEM 
elements  and  incorporate  a  special  assumed  stress  hybrid  EEM  formulation.  Incorporation  of  known  func¬ 
tional  forms  from  analytical  micromechanics  substantially  enhances  its  convergence.  A  high  level  of  accuracy 
with  significantly  reduced  degrees  of  freedom  has  been  achieved  with  VCEEM.  Computational  efficiency  is 
therefore  substantially  enhanced  compared  to  conventional  displacement-based  EE  models.  Successful  ap¬ 
plications  of  2D-small  deformation  VCEEM  have  been  made  in  thermo-elastic-plastic  problems  of  composite 
and  porous  materials  [61,  80].  An  adaptive  VCEEM  has  been  developed  in  [80],  where  optimal  improvement 
is  achieved  by  h-p  adaptation  of  the  displacement  field  and  p-enrichment  of  the  stress  field. 

The  cohesive  crack  propagation  model  has  been  incorporated  in  VCEEM  in  [35,  53]  to  model  interface 
debonding  in  fiber  reinforced  composites.  However,  in  these  models,  the  debonding  or  crack  evolution  path 
is  along  the  interface  and  hence  the  cohesive  zone  regions  are  known  a-priori.  In  the  event  that  the  crack 
branches  off  into  the  matrix,  the  path  is  no  longer  pre-assessed  and  needs  to  be  determined  at  each  load 
increment,  consistent  with  the  local  state  of  stresses,  strains  and  morphology.  This  task  is  considerably  more 
challenging  since  a  slight  deviation  can  lead  to  completely  wrong  prediction. 

The  motivation  of  this  work  is  derived  from  the  need  to  create  a  robust  finite  element  method,  extended 
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Voronoi  cell  finite  element  model  (X-VCFEM),  for  modeling  interface  debonding  with  arbitrary  crack  prop¬ 
agation  in  heterogeneous  materials.  This  is  an  essential  step,  prior  to  simulating  the  entire  microstructural 
failure  problem.  X-VCFEM  incorporates:  (a)  stress  discontinuities  across  the  cohesive  crack  through  branch 
functions  in  conjunction  with  level  set  methods,  (b)  crack  tip  stress  concentration  through  the  introduction 
of  multi-resolution  wavelet  functions  [38,  47,  71]  in  the  vicinity  of  the  crack  tip,  and  (c)  incremental  crack 
propagation  using  a  cohesive  energy  based  criterion  for  estimating  the  direction  and  length  of  the  incremental 
crack  advance. 

Finally,  a  Voronoi  cell  finite  element  model  is  developed  for  transient  elastodynamic  analysis  in  time 
domain  is  developed.  In  the  present  formulation,  the  inertia  field  is  approximated  in  terms  of  stresses  so 
as  to  satisfy  the  equilibrium  equation  a-priori.  The  weak  forms  of  kinematics  and  traction  reciprocity  are 
obtained  by  minimization  of  the  complementary  variational  principle.  Stress  wave  is  a  local  disturbance 
that  propagates  through  the  material,  resulting  in  high  stress  gradients  near  the  wave  front.  Therefore, 
localization  and  multi-resolution  properties  of  the  wavelet  functions  are  exploited  to  enhance  the  computa¬ 
tional  efficiency  by  enriching  the  stress  functions  only  locally  near  the  wave  front.  The  enrichment  is  carried 
out  adaptively  by  employing  posteriori  local  error  estimators  that  determine  the  required  translation  and 
dilation  of  the  wavelet  functions  at  each  time  step.  At  the  outset,  a  stable,  accurate  and  computationally 
efficient  adaptive  computational  framework  in  ID  is  developed  for  micro-mechanical  response  of  composites 
under  impact  loading.  The  accuracy  and  computational  efficiency  of  the  proposed  method  is  demonstrated 
through  comparison  with  analytical  solutions  and  conventional  FEM  packages. 

2.1  Organization  of  this  Report 

The  report  is  divided  into  five  subsequent  chapters.  In  chapter  3,  variational  formulation  and  various 
aspects  of  the  computational  scheme  for  stress  wave  propagation  in  composites  are  presented.  Extensions  of 
the  VCFEM  (X-VCFEM)  for  cohesive  crack  propagation  are  developed  in  Chapter  4.  Numerical  validation 
of  X-VCFEM  for  matrix  cracking  is  also  presented.  Based  on  the  preparation  of  previous  two  chapters. 
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the  X-VCFEM  for  modeling  interface  debonding  with  matrix  cohesive  cracking  are  developed  in  chapter  5. 
Finally,  an  account  of  multi-scale  modeling  is  presented  with  critical  examples. 
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Chapter  3 


Multi-resolution  Wavelet  Enriched 
Hybrid  Finite  Element  Method  for 
One  Dimensional  Elastic  Wave 
Propagation  in  Heterogeneous  Solids 

3.1  Introduction 

Heterogeneous  materials  are  being  used  increasingly  in  impact  related  applications  because  of  their  high 
strength  to  weight  ratio  and  improved  dynamic  properties  [32,  58].  As  a  consequence,  analysis  of  wave 
propagation  through  heterogeneous  media  is  being  pursued  consistently  by  researchers.  However,  most  of 
the  methods  available  in  the  literature  aim  at  obtaining  macroscopic  response  using  effective  properties  of 
the  composites.  The  major  challenge  in  obtaining  actual  micro-scopic  response  of  composites  under  dynamic 
loading  is  the  computational  size  of  the  problem  and  the  complex  nature  of  the  wave  propagation  phenomena. 
Stress  waves  experience  multiple  reflections,  transmissions  and  interference  in  the  presence  of  heterogeneities 
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in  the  microstructure  and  produce  dynamic  stress  concentrations  that  are  significantly  greater  than  that  in 
global  average  response  [69].  This  initiates  and  propagates  the  damage,  leading  to  failure  of  the  material. 
Recently,  attempts  have  been  made  to  account  for  the  micro-structural  effects  in  the  macro  response.  Wang 
and  Sun  [101]  developed  a  technique  that  includes  the  effect  of  micro-inertia  in  the  continuum  model  of  the 
heterogeneous  materials.  Using  these  macro-equations,  harmonic  and  transient  response  of  one-dimensional 
layered  medium  was  obtained  which  was  in  close  agreement  with  the  analytical  results  for  a  range  of  wave¬ 
lengths.  Fish  et  al  [46,  45]  developed  a  dispersive  model  for  wave  propagation  using  higher  order  homoge¬ 
nization  theory  with  multiple  spatial  and  temporal  scales.  A  goal-oriented  adaptive  modeling  technique  that 
solves  the  micro-mechanical  problem  using  actual  material  properties  was  developed  in  frequency  domain 
by  Romkes  and  Oden  [77].  The  method  used  local  error  estimators  to  identify  the  regions  where  actual 
properties  are  to  be  used  and  the  critical  frequencies  for  which  the  solution  is  to  be  improved. 

In  the  last  decade,  Voronoi  Cell  Finite  Element  Method  (VCFEM)  has  emerged  as  a  powerful  technique 
for  micro-mechanical  modeling  of  arbitrary  heterogeneous  materials.  VCEEM  has  been  developed  for  elastic, 
elasto-plastic  response  of  heterogeneous  materials  [61],  damage  initiation  and  propagation  in  ductile  as  well 
as  brittle  materials  [35,  36].  VCEEM  has  been  shown  to  be  significantly  more  efficient  than  the  conventional 
displacement  based  methods  for  2D  static  problems.  In  this  work,  a  hybrid  formulation  based  on  same 
principle  is  presented  for  one  dimensional  elastic  wave  propagation  in  heterogeneous  materials.  There  have 
been  a  few  hybrid/mixed  formulations  in  the  literature  for  elastodynamic  analysis  [3,  89,  31].  Inclusion  of 
adaptive  techniques  in  such  formulations  could  reduce  the  size  of  the  problem  substantially.  The  formulation 
proposed  in  this  work  utilizes  the  multi-resolution  properties  of  the  wavelets  for  adaptive  enrichment  of  stress 
function  so  that  least  number  of  stress  parameters  are  required. 

The  proposed  formulation  makes  two  independent  approximations:  stress  field  and  the  boundary  dis¬ 
placement  field.  The  internal  displacement  field  is  approximated  in  terms  of  stresses  so  as  to  satisfy  the 
equilibrium  equation  in  strong  sense.  The  kinematic  equation  in  the  element,  traction  reciprocity  and  the 
compatibility  of  the  internal  displacement  field  are  satisfied  in  weak  sense  as  Euler-Lagrange  equations. 
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The  chapter  is  organized  as  follows:  The  hybrid  stress  formulation  for  one  dimensional  wave  propagation 
is  introduced  first.  The  weak  form  is  derived  by  taking  variations  of  the  Hamiltonian  expressed  in  terms 
of  complimentary  energy.  The  stress  function  is  constructed  adaptively  using  multi-resolution  wavelets  and 
posteriori  error  indicator.  The  developed  adaptive  algorithm  is  implemented  for  simulating  wave  propagation 
through  layered  media. 

3.2  Assumed  Stress  Hybrid  Formulation 

A  typical  2D  Voronoi  cell  element  is  multi-phase  domain  consisting  of  inclusion/cavity  surrounded  by  matrix 
phase  [61].  However,  the  two-noded  ID  element  presented  here  consists  of  a  single  phase  such  that  each  layer 
of  heterogeneous  layered  media  becomes  one  element  fie  without  need  of  any  further  refinement.  The  element 
boundary  9fle  with  unit  outward  normal  n®  is  formed  by  two  nodes  and  may  consist  of  prescribed  traction 
Tie,  prescribed  displacement  r„e  and  inter-element  boundary  Tme-  For  the  hybrid  element  formulation,  in 
the  absence  of  body  forces,  the  micro  mechanics  elastodynamic  initial  boundary  value  problem  is  described 
as: 


Find  ((T,  u,  Ur)  G  T  x  V  x  Vr  satisfying 


dB 

{a) 

V  •  (T  =  pii  and 

^  =  e  G  De 
da 

Ur  =  u  on  Tue  , 

CT  •  n®  =  t  on  Fie 

(b) 

(3.1) 

<7  =  (To  ,  u  =  Uo 

and  ii  =  iio  at  t  =  0 

(c) 

where  V  =  The  variables  cr,  e  and  B  are  the  equilibrated  stress  field,  the  corresponding  strain  field  and 
the  complimentary  energy  respectively  in  the  element  interior.  T,  V  and  Vr  correspond  to  Hilbert  spaces 
containing  the  stress,  internal  displacement  and  boundary  displacement  solutions  respectively.  ur{x,t)  is 
kinematically  admissible  compatible  displacement  field  at  the  element  boundary  (nodal  displacement),  while, 
u{x,  t),  u{x,  t)  and  u{x,  t)  represent  kinematically  admissible  internal  displacement,  velocity  and  acceleration 
fields  respectively.  Furthermore,  p  denote  the  material  density  and  i{t)  is  the  applied  traction  (F(t)  =  i{t)A 
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is  nodal  force  where  A  is  cross-sectional  area). 
The  Hamiltonian  for  one  element  is  written  as 


ff(t)  =  [  '  (T(t)  -  V{t))dt  (3.2) 

Jti 

where  T{t)  is  the  Kinetic  Energy  and  V{t)  is  the  Potential  Energy.  Since  the  internal  displacement  field 
u  is  not  compatible,  the  incompatible  displacement  at  the  boundary  can  be  made  zero  by  including  the 
condition  tt  —  itr  =  0  as  a  constraint  to  the  Hamiltonian  and  applying  the  Lagrange  Multiplier  technique. 
Also,  expressing  the  Potential  Energy  in  terms  of  Complementary  Energy  density,  and  applying  divergence 
theorem,  the  Hamiltonian  can  be  rewritten  as 


■  dt 


H{t)  =  f  /  piiud^l  +  [  Bd^+  f  VCTitdflli 
Jti  Jfi  Jfi  Jfi  J 

+  f  f  anudd^l  +  f  turdd^l  +  f  t{u  —  UT)dd^l\ dt 

Jti  I  Jan^  Jan^  Jan^  J 


(3.3) 


where,  t  =  an  is  boundary  traction.  In  the  hybrid  formulation,  the  equilibrium  and  constitutive  equations, 
and  the  displacement  compatibility  at  the  inter-element  boundary  (i.e.  nodes)  are  satisfied  a-priori  in 
strong  sense.  The  kinematic  equation,  traction  boundary  condition,  traction  reciprocity  on  the  inter-element 
boundary  and  compatibility  of  internal  displacement  field  are  obtained  as  the  Euler-Lagrange  equations  from 
the  stationarity  of  the  complimentary  energy  functional  (3.3)  as  follows:  Setting  the  first  variation  of  energy 
functional  (3.3)  equal  to  zero 


SH  =  f  f  piiSud^l  +  [  dadfl -I-  [  V{Sa)ud^l+  [  VaSud^lXdt  (3.4) 

Jti  I  Jq.  Jq.  oa  Jq  Jq  j 

+  j  {—I  SanuddQ,  +  I  tSurddU  +  j  5an{u  —  itr)d9fl  —  /  anSurddU  >  dt  =  0 

Jti  I  Jan^  Jan^  Jan^  Jan^  ) 


where 


\^It  ~  ^  I 


has  been  used  [102].  Since  equation  (3.5)  is  valid  for  any  arbitrary  time  interval  {ti,t2},  the  term  inside  the 
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time  integral  can  be  set  to  zero.  Applying  divergence  theorem  to  the  sixth  term  in  equation  (3.5),  and  noting 
that  equilibrium  and  constitutive  equations  are  satisfied  a-priori,  we  arrive  at  the  following  weak  form: 


/  (e  —  Vit)  SadU  +  I  it  —  an)  SurddU  +  I  San  (it  —  itr) 

Jn  Jan  Jan 


ddU  =  0 


(3.5) 


which  results  in  weak  satisfaction  of  following  equations 


e  =  Vit 

in 

Oe 

(Kinematics) 

an  =  t 

on 

dUt 

(Traction  Boundary) 

an'^  = 

an~ 

on  d^lm 

(Traction  Reciprocity) 

u  =  Ur 

on 

dUg 

(Compatibility  of  internal  displacements) 

(3.6) 

3.3  Hybrid  Element  Assumptions  and  Weak  Form 

3.3.1  Element  Assumptions 

In  the  hybrid  stress  formulation,  stress  field  in  the  domain  is  approximated  as 

a{x,  t)  =  [P(x,  t)]  {/3{t)}  (3.7) 

where  [P(x,t)]  is  the  matrix  containing  stress  interpolation  functions  that  are  functions  of  time,  and 
are  the  corresponding  unknown  coefficients.  To  minimize  the  computational  cost,  it  is  desirable  to  have 
least  number  of  terms  in  the  stress  interpolation  function.  In  the  present  formulation  the  inertia  field  in  the 
domain  is  approximated  in  terms  of  stress  approximation.  The  acceleration  field  is  interpolated  in  such  a 
way  that  the  equilibrium  equation  is  satisfied  a-priori. 

u{x,  t)  =  ^  t)]  (3-8) 
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The  velocity  and  internal  displacement  fields  in  the  domain  are  obtained  by  numerical  integration  of  accel¬ 
eration.  According  to  Newmark  method,  which  is  an  implicit  scheme,  velocity  and  displacement  are  given 
as 


^t+At  ^ 

asu*  +  a4U*  +  051**“*“^* 


(3.9) 

(3.10) 


where  cti  =  (1  —  6)At,  a2  =  SAt,  =  At,  =  (1/2  —  a)At^  and  =  aAt^,  where  a  and  6  are  the  inte¬ 
gration  constants  and  superscript  ()*  denotes  quantities  at  previous  time  step.  The  boundary  displacements 
are  approximated  independently  as 

{itr(x,  t)}  =  [L(x)]  {qr(f)}  (3.11) 

In  this  one  dimensional  formulation,  the  boundary  displacements  {itr}  are  same  as  nodal  displacements 
{qr}  and  therefore  [L]  is  an  identity  matrix. 

3.3.2  Weak  Form 

Substituting  equilibrium  and  constitutive  equations  in  equation  (3.5),  the  weak  form  of  complimentary  energy 
functional  can  be  written  as 

I  6aSad^l+  j  V{6a)ud^l+  j  t6urdd^l+  j  danurdd^l  —  j  an6urdd^l  =  0  (3.12) 

Jn  Jn  Jan^  Jan^  Jan^ 

The  matrix  equation  of  the  weak  form  (3.12)  is  obtained  by  substituting  the  approximations  for  stress, 
internal  displacement  and  boundary  displacement  fields  from  equations  (3.7),  (3.10)  and  (3.11)  into  equation 
(3.12) 

{Sfdf  [H]  {/3}  +  {Sfdf  {i^4  +  as  {Sfdf  [M]  {/3} 

+  {Rtf{Sqr}  -  {Sfdf[G]{qr}  -  {fdf[G]{Sqr}  =  0  (3.13) 
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where  [H]  =  [P]^[S][P]dO  [M]  =  i  [VP]^[VP]ciO 

[G]  =  Ian.  [P]^[n][L]ci50  {R*}  =  {t}^[n][L]ciri„ 

{Rj}  =  [VP]^  {{it}*  +  az{uy  +  aiiuY}  dU  (3.14) 

3.4  Construction  of  Wavelet  Based  Stress  Functions 

Choice  of  functions  for  stress  interpolation  is  the  most  imperative  task  in  the  hybrid  stress  formulations  as 
it  determines  the  computational  effort  required  and  the  level  of  accuracy  that  can  be  achieved.  VCFEM 
formulations  in  the  past  [61,  54]  have  employed  polynomial,  reciprocal  and  some  special  functions  like  branch 
functions,  wavelets  etc  for  stress  interpolation  to  enhance  the  accuracy  and  computational  efficiency. 
Transient  wave  propagation  in  elastic  solids  is  essentially  traverse  of  a  local  disturbance  of  high  stress 
gradients  through  the  material.  The  formulation  was  tested  with  stress  interpolation  based  on  polynomial 
functions.  It  is  observed  that  polynomials  of  order  as  high  as  15  are  not  able  to  capture  the  abrupt  variations 
in  the  stress  field  at  the  moving  wave  front.  Wavelets  are  the  functions  that  have  localization  and  multi¬ 
resolution  properties  which,  when  coupled,  facilitate  local  enrichment  of  the  stress  function  and  therefore  are 
the  most  suitable  candidates  for  this  purpose.  A  brief  introduction  to  wavelets  is  provided  in  the  following 
paragraphs,  which  is  followed  by  details  on  construction  and  use  of  wavelet  basis  for  stress  wave  propagation 
problems. 

3.4.1  Principles  of  Wavelets  and  Multi-resolution  Analysis 

Scaling  function  and  wavelet  function  ^(x)  are  the  basic  building  blocks  of  the  multi-resolution  analysis. 
Scaling  function  is  defined  as  a  recursive  function  that  satisfies  the  two-scale  relation 

^(x)  =  ^p(fc)^(2x  —  k)  (3.15) 

k 

where  {p{k)}kez  are  the  filter  coefficients.  The  scaling  function  has  a  compact  support  if  only  a  finite  number 
of  coefficients  p{k)  are  non-zero.  Translation  of  scaling  function  ^  by  a  factor  of  2"  and  dilation  by  a  factor 
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of  fc  •  2  "  forms  unconditional  basis  of  subspace  y„  C  (TZ)  as 

(x)  =2"/2^(2"x-fc)  (3.16) 

where  n  is  resolution  level.  The  scaling  function  (j)  is  orthonormal  if  translations  at  the  same  resolution  level 
satisfy  the  orthogonality  condition 


/: 


<pn,k{x)'Pn,i{x)dx  =  6k, I  V  n,k,l  €  Z 


(3.17) 


If  the  scaling  function  is  orthonormal,  the  best  approximation  of  a  function  f{x)  at  resolution  level  n  is 
expressed  as  the  orthogonal  projection  of  /  on  subspace  Vn  as: 


/OO 

f{x)(j)n,k{x)dx  (3.18) 

K 

In  general,  approximation  of  /(x),  at  resolution  level  n  is  contained  in  approximation  at  any  resolution  level 
higher  than  n  i.e.  {0}  =  y_oo  C  •  •  •  C  y-i  C  Vb  C  y  C  •  •  •  C  yoo  =  LP‘{TZ).  This  means  that  function 
/(x)  is  approximated  better  at  higher  resolutions  and  some  information  is  lost  in  transition  from  higher  level 
Vn+i  to  lower  level  y,.  This  difference  is  characterized  by  an  orthogonal  complementary  subspace  Wn  so 
that  Vn+i  =  y  ©  Wn  V  n.  A  basis  that  spans  the  subspace  Wn  can  be  obtained  in  the  same  manner  as  for 
scaling  function,  i.e.  by  translation  and  dilation  of  the  mother  wavelet  function 


^/)(x)  =  q{k)'tp{2x  -  k) 
k 


(3.19) 


The  wavelet  basis  is  orthonormal  if  any  two  translated  and/or  dilated  wavelets  satisfy  the  orthogonality 
condition 


/: 


i>n,k{x)'tpni,i{x)dx  =  5n,ni6k,i  V  n,m,k,l£  Z 


(3.20) 
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The  wavelet  basis  is  semi-orthogonal  if  any  two  translated  wavelets  at  different  resolution  levels  satisfy  the 
semi-orthogonality  condition 


/: 


'>pn,k{x)'>pm,i{x)dx  =  0,  Ti  =  TTi  V  n,m,k,l  £  Z 


(3.21) 


An  approximation  of  the  function  f{x)  at  the  n  —  th  resolution  level  may  be  expressed  as  the  orthogonal 
projection  of  /  on  Wn  as 


/OO 

f{x)'ll)n,k{x)d3. 

-OO 


(3.22) 


Thus,  approximation  of  function  f{x)  at  higher  resolution  can  be  obtained  as 


'^n+lf{x')  —  ^  ^^(in,k^n,k{x')  -|-  ^  ^  (^)  (3.23) 

k  k 

These  multi-resolution  properties  of  wavelet  functions  provide  the  basis  for  adaptive  enrichment  in  the 
regions  where  residual  is  higher  at  the  lower  resolution  level. 

3.4.2  Selection  of  the  Wavelet  Function 

In  the  present  formulation,  approximation  for  acceleration  field  is  constructed  on  derivative  of  the  stress 
interpolation  functions.  Also,  the  calculation  of  error  norm  discussed  in  section  (3.4.4)  involves  second 
derivatives  of  stress  functions.  Therefore,  it  is  desirable  that  the  stress  functions  be  differentiable  and  have 
explicit  analytical  expressions.  One  of  the  most  commonly  used  wavelet  functions  is  Daubechies’  compactly 
supported  orthonormal  wavelets  [21,  23].  However,  they  are  constructed  through  recursive  algorithms  and 
do  not  have  explicit  analytic  expressions,  therefore  are  not  suitable  for  present  formulation.  On  the  other 
hand,  Chui- Wang’s  B-spline  wavelet  bases  [16,  62]  have  explicit  analytic  expressions  for  scaling  and  wavelet 
functions,  and  therefore  are  implemented  in  this  formulation.  Chui- Wang  wavelets,  which  are  semi-orthogonal 
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and  compactly  supported,  are  based  on  B-spline  functions,  which  are  defined  by  recursive  convolution 

/oo  rx 

Nto_i(x  -  t)Ni(t)dt  =  /  (3.24) 

-OO  ^  x  —  1 

where  Ni(x)  is  a  box  function.  The  two  scale  relation  for  this  scaling  function  is  given  as 

m  I 

cf>{x)  =  N„(x)  =  -  k)  (3.25) 

The  corresponding  wavelet  basis,  which  satisfies  the  semi-orthogonality  condition  (3.21),  is  given  by 

2m— 2 

v>(x)  =  2-™+i  Y.  (-l)''N2™(fc  +  1)n(™)(2x  -  k)  (3.26) 

fc=0 

where 

m  I 

=  (3.27) 

For  B-spline  wavelet  bases,  the  scaling  function  and  the  wavelet  function  are  compactly  supported  i.e.  they 
are  defined  on  a  finite  closed  interval.  The  support  for  scaling  and  wavelet  functions  are  given  by 

supp  =  [2^"fc,2^"(fc -l-m)] 

supp  =  [2^"fc,  2^"(fc -I- 2m  —  1)]  (3.28) 

The  finite  number  of  non-zero  translations  that  form  a  basis  for  interpolation  on  an  interval  [a,  b]  can  be 
identified  using  the  above  expressions  for  the  support  as 

(2"a  -  m  -I- 1)  <  <  (2"6  -  1)  (3.29) 

(2"a-m-l-l)  <  <  (2"6- 1)  (3.30) 

The  total  non-zero  translations  over  the  interval  [a,  b]  for  the  scaling  function  and  the  wavelet  function  are 
[2"  (6  -  o)  -I-  m  -  1]  and  [2"  (6  -  a)  -I-  m  -  1]  respectively. 
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3.4.3  Multi-resolution  Wavelet  Based  Stress  Functions 


The  wavelet  based  stress  function  is  constructed  by  forming  a  basis  by  translating  the  scaling  function  at 
resolution  n.  The  required  translations  to  span  the  volume  of  the  element  are  calculated  using  expression 
(3.29).  The  stress  function  is  expressed  as 


CT  =  ^  4)n,kix)^kit)  (3.31) 

The  stress  function  is  enriched  locally  in  the  vicinity  of  wave  front  by  using  wavelet  functions  at  resolution 
n  to  form  a  basis  at  increased  resolution  level  (n  +  1) 


CT  =  '^4>n,kix)Pkit)  +  '^‘il’n,i{x)Pi{t)  (3.32) 

k  I 

The  local  region  can  be  enriched  further  to  resolution  level  (n  +  2)  by  using  wavelet  functions  at  resolution 
(n  +  1)  and  so  on  until  the  desired  accuracy  is  achieved.  While  forming  the  wavelet  basis  in  the  element, 
some  translations  of  the  scaling  and  wavelet  functions  fall  partially  outside  the  element.  In  such  cases,  the 
wavelets  are  truncated  at  the  element  boundary. 


3.4.4  Error  Criteria  for  Time  Dependent  Adaptive  Enrichment 


It  is  imperative  to  accurately  determine  the  location  of  the  wave  front  where  the  local  enrichment  of  stress 
function  is  to  be  carried  out,  for,  it  governs  the  accuracy  of  the  solution  and  also  determines  the  computational 
cost.  As  seen  in  section  (3.2),  the  kinematic  equation  (3.6)  in  the  domain  is  satisfied  in  an  average  sense. 
Therefore,  the  residual  in  satisfying  this  equation  is  used  as  error  indicator  for  adaptive  enrichment.  The 
error  norm  is  defined  as 


1 

AO  leL  _ 


X  100 


(3.33) 


where  AO  and  \e\max  are  volume  of  subdomain  and  maximum  strain  in  the  element  respectively.  The  domain 
is  discretized  into  2"  number  of  subdomains  and  error  norm  (3.33)  is  calculated  in  each  subdomain.  The 
subdomains  for  which  the  norm  exceeds  a  pre-defined  tolerance,  are  the  regions  where  enrichment  is  to  be 
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carried  out  using  the  procedure  outlined  in  previous  section.  This  increases  the  resolution  level  to  {n  +  1). 
The  posteriori  error  analysis  and  corresponding  enrichment  is  repeated  until  the  tolerance  is  met. 


3.5  Solution  Method  and  Numerical  Aspects 

3.5.1  Solution  for  the  Field  Variables 

The  matrix  equation  (3.13)  of  the  weak  form  of  complimentary  energy  principle  can  be  rewritten  as 

{[H]{/3}  +  {Ri}  +  a5[M]{/3}  -  [G]{qr}}  +  {{Rtf  -  {/3}^[G]}  {Sqr}  =  0  (3.34) 

Since  {6/3}  and  {dqr}  are  arbitrary  and  can  be  varied  independently,  and  the  corresponding  bracketed  terms 
are  independent  of  these  variations,  the  two  bracketed  terms  should  vanish  individually.  Setting  the  first 
bracketed  term  in  equation  (3.34)  equal  to  zero  gives  the  local  equations  for  each  element: 

[HM]e{/3}e  =  [G]e{gr}e  "  {R/}e  (3.35) 

where 

[HmIb  =  [H]e  -  a5[M]e  (3.36) 

Adding  energy  of  all  elements,  and  setting  the  second  bracketed  terms  equal  to  zero  gives  weak  form  of  global 
traction  reciprocity  condition 

=  (3.37) 

e=l  e=l 

If  the  element  [HmIb  matrix  is  invertible,  the  stress  coefficients  can  be  expressed  in  terms  of  nodal  displace¬ 
ments  using  equation  (3.35).  The  static  condensation  of  equations  (3.35)  and  (3.37)  gives  linear  system  of 
simultaneous  equations 

E[G]nHM]-MG]e{gr}e  =  (3-38) 

e=l  e=l 

or  [K]{qr}  =  {F}  (3.39) 
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which  can  be  solved  for  the  nodal  displacement  vector  {qi?}. 


3.5.2  Stability  Conditions 

The  stability  conditions  for  multi-field  mixed  variational  formulations  are  derived  in  [4,  12,  104].  Within 
this  framework,  the  stability  conditions  for  the  stress-displacement  field  variational  problem  in  the  dynamic 
hybrid  FEM  can  be  stated  as  follows: 

•  The  matrix  [Hm]  should  be  positive  definite.  This  also  ensures  invertability  of  the  [Hm]  matrix.  From 
the  definition  of  [Hm]  in  equation  (3.35),  the  necessary  condition  for  it  to  be  positive  definite  is 

{xY[liM]{x}  >0  {x}^[H]{x}  >  0  &  {x}^[M]{x}  >  0 

i.e.  the  matrix  [H]  should  be  positive  definite  and  matrix  [M]  should  be  positive  semi-definite.  For 
matrix  [H]  to  be  positive  definite,  firstly,  [S]  be  positive  definite,  which  is  true  for  elastic  problems. 
Secondly,  the  finite-dimensional  stress  subspace  T  should  be  spanned  uniquely  by  the  basis  functions 
[P].  This  is  satisfied  by  assuring  linear  independence  of  the  columns  of  basis  functions  [P].  As 
discussed  in  previous  section,  stress  function  contains  a  basis  formed  by  low  resolution  scaling  function 
everywhere  in  the  domain,  and  is  enriched  locally  using  wavelet  function.  Though  the  wavelet  functions 
are  orthogonal  to  the  scaling  function,  the  orthogonality  is  destroyed  near  the  element  boundaries  where 
scaling  and  wavelet  functions  are  partially  outside  the  element.  In  situations  where  only  a  small  portion 
of  the  scaling  and  wavelet  functions  fall  inside  the  element,  these  portions  could  become  dependent  or 
nearly  dependent.  To  encounter  this  problem,  the  rank  of  the  [P]  matrix  is  first  determined  from  the 
diagonal  matrix  resulting  from  a  Cholesky  factorization  of  the  square  matrix 

[H*]  =  [  [P]^[P]dO 

Nearly  dependent  columns  of  [P]  will  result  in  very  small  pivots  during  Cholesky  factorization.  The 
corresponding  wavelet  function  terms  are  dropped  from  the  stress  function  to  prevent  numerical  inac- 
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curacies  in  inverting  [Hm]- 


•  To  ensure  non-zero  stress  field  in  the  element  for  all  non-rigid  body  displacement  fields  on  the  element 
boundary  uf ,  the  dimensions  of  the  stress  and  displacement  subspaces  must  satisfy  stability  condition 
rip  >  rig^  —  1,  where  np  is  the  number  of  /3  parameters,  and  rig^  is  the  number  of  displacement  degrees 
of  freedom  on  the  element  boundary.  In  this  one  dimensional  formulation,  this  is  always  satisfied  if 
np  >  1. 

Another  important  factor  that  determines  the  stability  and  accuracy  of  solution  of  dynamic  problems  is 
the  stability  of  the  time  integration  scheme.  Newmark  method,  being  used  here,  is  unconditionally  stable  and 
most  accurate  when  the  integration  constants  in  equations  (3.10)  and  (3.9)  take  values  d  =  |  a  =  j,  which 
is  constant-average-acceleration  method,  also  called  trapezoidal  rule  [6].  Therefore,  there  is  no  minimum 
time  step  size  requirement  and  it  is  the  accuracy  of  the  solution  that  decides  the  time  step  size.  In  wave 
propagation  problems,  the  time  step  size  is  determined  as  [6] 

At=^ 

nc 

where  is  the  critical  wavelength  to  be  represented,  n  is  the  number  of  time  steps  necessary  to  represent 
the  travel  of  the  wave,  and  c  =  is  the  wave  speed. 

3.6  Numerical  Examples 

3.6.1  Wave  Propagation  through  Layered  Media 

The  problem  of  wave  dispersion  in  layered  media  presented  in  [46]  is  considered  here.  A  layered  bar  composed 
of  two  materials  with  properties  Ei  =  200  GPa,  E2  =  5  GPa  and  pi  =  P2  =  8000kg/m^  is  shown  in  figure 
3.1.  The  lengths  of  the  two  layers  in  a  unit  cell  are  h  =  I2  =  0.01  m  and  there  are  50  unit  cells.  One  end 
of  the  bar  is  fixed  while  other  end  is  subjected  to  impact  load  F{t)  =  Fq  —  T)^[l  —  h{t  —  T)], 

where  T  =  15.71  ps  is  the  duration  and  Fq  =  50  KN  is  the  amplitude  of  the  impact  load,  and  h{t)  is 
heaviside  function.  Response  of  ABAQUS  model  with  50  T2D2  elements  per  layer  is  taken  as  a  reference 
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and  the  response  of  hybrid  finite  element  model  with  one  element  per  layer  is  compared.  Figure  3.4  shows 
displacement  at  x  =  0.5  m  as  a  function  of  time.  It  can  be  observed  that  proposed  formulation  predicts  the 
phenomenon  of  dispersion  very  accurately. 

3.6.2  Effect  of  Assumption  of  Periodicity 

The  homogenization  methods  often  make  assumption  of  periodicity  in  order  to  include  microstrucutral  effects 
in  the  macro-response.  In  this  example,  the  accuracy  of  such  assumption  is  investigated.  A  composite  bar 
with  properties  given  in  previous  example  and  with  length  li  =  I2  =  0.1  m  is  considered.  Figure  3.2  (a)  and 
(b)  show  a  unit  cell  with  periodic  boundary  conditions  and  a  full  model  of  50  layers  with  fixed  end  conditions 
respectively.  In  this  case  T  =  63  jis  is  the  duration  and  Fq  =  100  KN  is  the  amplitude  of  the  impact  load. 
Figures  3.5  and  3.6  show  stress  response  near  the  boundary  and  away  from  the  boundary  as  a  function  of  time 
respectively.  It  can  be  observed  that  interference  of  waves  reflected  from  the  boundary  produce  significantly 
different  response  near  the  boundary.  The  assumption  of  periodicity  is  reasonably  accurate  away  from  the 
boundary.  However,  response  deviates  as  the  reflected  waves  arrive  and  interfere  with  the  incoming  waves. 

3.7  Conclusions 

An  assumed  stress  hybrid  Voronoi  cell  finite  element  model  for  analysis  of  elastic  wave  propagation  in  het¬ 
erogeneous  materials  is  proposed  in  this  work.  Stress  field  in  the  domain  and  compatible  displacement 
field  at  the  boundary  are  interpolated  independently.  The  nonconforming  internal  displacement  field  is  ap¬ 
proximated  in  terms  of  stresses  such  that  equilibrium  equation  is  satisfied  pointwise.  Stress  functions  are 
based  on  low  resolution  B-spline  scaling  functions  which  are  adaptively  enriched  using  wavelet  functions  in 
the  local  region  of  high  stress  gradients  determined  by  residual  based  posteriori  error  indicator.  Adaptive 
enrichment  of  the  stress  function  exploiting  multi-resolution  properties  of  wavelet  functions  reduces  the  de¬ 
grees  of  freedom  of  the  problem  and  enhances  the  computational  efficiency  significantly.  As  demonstrated 
through  comparison  with  standard  FEM  packages,  the  proposed  formulation  predicts  the  phenomenon  of 
wave  reflection,  transmission,  dispersion  etc  accurately. 


29 


This  work  advocates  multi-resolution  wavelet  enriched  hybrid  FEM  as  a  potential  method  for  micro-mechanical 


response  of  composites  under  impact  loading  in  one  dimension.  Wave  propagation  through  heterogeneous 
solids  in  two  dimensions  is  more  involving  due  to  complex  interaction  of  dilatational  and  distortional  waves 
at  the  boundaries  of  heterogeneities.  In  the  subsequent  work,  this  formulation  will  be  extended  in  two 
dimensions  to  develop  dynamic  VCFEM  for  investigating  propagation  of  elastic  waves  in  heterogeneous 
microstructures. 
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Figure  3.1:  Wave  propagation  through  layered  media:  (a)  Unit  cell,  (b)  Impact  load  as  a  function  of  time, 
(c)  Composite  bar 
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Figure  3.2:  Effect  of  assumption  of  periodicity:  (a)  Periodic  model,  (b)  Full  model 
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Figure  3.3:  Displacement  response  at  the  center  of  a  layered  bar:  Comparison  with  Abaqus 


Figure  3.4:  Displacement  response  at  the  center  of  a  layered  bar:  Comparison  with  Abaqus 


Chapter  4 


The  Extended  Voronoi  Cell  Finite 
Element  Model  For  Multiple  Cohesive 
Cracks  Propagation 

4.1  Introduction 

Numerical  analysis  and  simulation  of  the  growth  of  multiple  cracks  in  materials  is  a  challenging  enterprise 
due  to  morphological  and  constitutive  complexities  that  govern  its  growth.  The  conventional  finite  element 
method  suffers  from  very  slow  convergence  since  the  element  formulation  does  not  account  for  high  gradients 
and  singularities.  Even  a  very  high  density  mesh  cannot  overcome  pathological  mesh  dependence  near  the 
crack  tips  and  avoid  biasing  the  direction  of  crack  propagation. 

In  this  chapter,  an  extended  VCFEM  or  X-VCEEM  is  developed  for  modeling  the  growth  of  multiple  cohe¬ 
sive  cracks  in  a  brittle  material.  The  model  accounts  for  interaction  between  cracks  and  invokes  an  adaptive 
crack  growth  formulation  to  represent  the  continuously  changing  direction  of  evolving  cracks.  X-VCEEM 
augments  the  conventional  VCEEM  model  by  incorporating  multi-resolution  wavelet  functions  [38,  47,  71] 
in  the  vicinity  of  the  crack  tip,  in  addition  to  branch  functions  based  on  level  set  methods.  The  incremental 
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crack  propagation  direction  and  length  are  adaptively  determined  by  a  cohesive  energy  based  criterion.  No 
remeshing  is  needed  in  X-VCFEM  for  simulating  crack  growth,  and  this  adds  to  its  desirability  and  effec¬ 
tiveness.  It  begins  with  the  X-VCFEM  formulation,  followed  by  numerical  example  showing  the  convergence 
of  this  model.  Then,  X-VCFEM  is  used  to  understand  the  influence  of  cohesive  parameters,  e.g.  peak 
stress  and  critical  separation  on  crack  growth  in  a  monolithic  brittle  material.  Subsequently,  the  effect  of 
morphological  distributions  including  crack  interaction,  clustering,  alignment,  etc.  on  growth  and  merging 
are  studied  as  important  factors  critical  to  the  failure  process. 

4.2  Voronoi  Cell  Fern  Formulation  for  Multiple  Propagating  Cracks 

The  Voronoi  cell  finite  element  mesh  for  a  brittle  matrix  with  a  dispersion  of  pre-existing  cracks  is  shown 
in  figure  4.1(a).  The  typical  Voronoi  cell  mesh  corresponds  to  an  unstructured  mesh  that  is  generated  by 
Dirichlet  or  Voronoi  tessellation  of  the  domain,  based  on  the  position,  shape  and  size  of  heterogeneities  (in¬ 
clusion,  void,  crack  etc.).  Various  tessellation  schemes  have  been  discussed  and  developed  in  [35,  61].  While 
the  name  Voronoi  cell  has  been  historically  used  because  of  its  association  with  point  seeds  in  the  generation 
process,  the  cells  used  in  VCFEM  may  be  variants  of  this  construct.  Essentially  they  represent  neighbor¬ 
hood  or  regions  of  influence  for  each  heterogeneity.  Subsequently  the  Voronoi  cell  EE  formulation  considers 
each  cell  as  a  super-element  consisting  of  a  heterogeneity  and  its  neighborhood  surrounding  matrix  [61,  80] 
without  any  further  subdivision.  The  interfacial  debonding  analyses  in  [35,  53]  invoke  the  cohesive  zone 
models  to  represent  the  growth  of  interfacial  crack.  However  the  main  difference  between  that  formulation 
and  the  present  one  is  that,  in  the  present  case  the  path  of  the  crack  is  arbitrary  and  is  a-priori  unknown. 
This  poses  significant  challenges  that  have  been  overcome  with  the  X-VCFEM  formulation. 

Consider  a  pre-cracked  microstructural  region  fl  consisting  of  N  cracks  as  shown  in  figure  4.1(a).  The 
region  is  divided  into  an  unstructured  finite  element  mesh  of  arbitrary  Voronoi  cells.  A  typical  VC  element 
fie  containing  a  crack  and  its  neighboring  matrix  is  depicted  in  figure  4.1(b).  The  element  boundary  9fle  with 
outward  normal  may  consist  of  regions  of  prescribed  traction  Eje,  prescribed  displacement  E^e  and  inter- 
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element  edges  i-e.  =  rieUr^eUrme-  Furthermore,  each  element  consists  of  a  crack  containing 
a  fracture  process  zone  that  is  represented  by  a  cohesive  zone  model.  The  incompatible  displacement  field 
across  the  crack  Ter  is  facilitated  through  a  set  of  connected  node-pairs  along  the  crack  length.  The  node-pair 
merges  at  the  crack  tip  by  enforcing  the  same  displacement.  The  normal  along  the  crack  path  is  denoted  by 
ncr  Pqj.  VCFEM  element  formulation,  the  micromechanics  boundary  value  problem  is  described  as: 

Find  ((T,u^,u"’’)  G  r  X  X  V"’’  satisfying 

-  dB 

V  •  (7  -I-  f  =  0  and  =  e  G  fie  (a) 

0(7 

u"®  =  u  on  F„e  ,  (7  •  n"®  =  t  on  Fje  and  cr  •  =  1^°'^  on  T^r  {b)  (4.1) 


The  variables  (7,  e,  B  and  f  are  the  equilibrated  stress  fields,  the  corresponding  strain  fields,  the  complimen¬ 
tary  energy  and  body  forces  per  unit  volume  respectively  in  the  element  interior.  T,  and  correspond 
to  Hilbert  spaces  containing  the  stress  and  displacement  solutions  respectively,  is  the  kinematically 
admissible  displacement  field  on  the  element  boundary  dTtf  and  represents  the  displacements  on  the 
internal  cohesive-crack  surfaces  Ter-  Variables  with  superscript  E  are  on  the  element  boundary  while  those 
with  superscripts  cr  correspond  to  the  crack  surface.  The  traction  between  node-pairs  on  the  crack 
surface  are  modeled  by  the  cohesive  zone  traction-separation  law.  The  VCFEM  formulation  is  based  on 
the  assumed  stress  hybrid  finite  element  method,  in  which  stationarity  conditions  of  the  element  energy 
functional  in  the  variational  principle  yields  weak  forms  of  the  kinematic  equation  and  traction  reciprocity 
conditions,  as  Euler  equations.  In  the  small  deformation  elasticity  incremental  formulation  for  evolving 
cracks,  the  element  energy  functional  He  is  defined  in  terms  of  increments  of  stresses  and  displacements  as: 


-b 

-b 


lle{aij,Aaij,uf,Auf,Ui'^,Aui'^)  =  -  AB{aij,Aaij)d^l  -  /  eijAaijdi 
[  {aij  +  Aaij)nf{uf  +  Auf)ddTl  -  f  {U  +  Ati){u^i  +  Au^i)dVte 

JVtm 

(uy  +  Aai^)nf{uf  +  Auf)dTer  -  j,  +  Aai^)nf{uf  +  Auf)dT. 

IJl 


-I- A-u?'' A-ur  1  2 

tr’^diu^r  -uT)dTe 


(4.2) 
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where  B  =  |(t  :  S  :  cr  is  the  complimentary  energy  density  and  AB{aij,  Aaij)  is  its  increment  due  to  stress 

1  2 

increase.  S  is  the  material  compliance  matrix.  The  notations  (•)  and  (•)  represent  two  sides  of  the  internal 
cohesive  crack  surface.  The  last  term  provides  the  work  done  by  the  cohesive  tractions  due  to  crack 
surface  separation.  In  VCFE  formulation,  the  equilibrium  conditions  and  constitutive  relations  in  the  matrix 
and  the  compatibility  conditions  on  the  element  boundary  and  crack  surface  are  satisfied  a-priori  in  a  strong 
sense.  The  element  kinematic  equation: 

VUe  =  ee  in  fie  (4.3) 

is  however  satisfied  in  a  weak  sense  from  the  stationary  condition  of  the  element  energy  functional  in 
equation(4.2).  The  weak  form  is  obtained  by  setting  the  first  variation  of  He  with  respect  to  stress  in¬ 
crements  to  zero,  i.e. 


jfle  dAdij 


-feij)  SAdij  dfl  -l- 


SAaij  rij  (uf  +  Auf)  ddUg 


+  j ^  6 Aaij  nf  (uf  +  Au^)  dV„  -  j,  SAaij  nj  (4"  +  Au^)  dV,r  =  0 


(4.4) 


Solution  of  equation  (4.4)  yields  domain  stresses.  Furthermore,  the  VCFE  formulation  assumes  weak  sat¬ 
isfaction  of  the  traction  reciprocity  conditions  on  (i)  the  inter-element  boundary  T^e,  and  (iii)  the  domain 

1  2 

traction  boundary  Fje  and  (iii)  the  crack  surfaces  F^  and  F^: 


[aij  +  Aaij)n^'^  =  —{aij+Aaij)n^^  on  F^e  (inter-element  boundary) 
{aij  +  Aaij)n^  =  ii+Aii  on  F^e  (traction  boundary) 

{oij  +  AaijYnJ  =  {oij  +  AaijYnJ  on  Fcr 


(4.5) 


In  the  variational  principle,  the  weak  form  is  obtained  by  setting  the  first  variation  of  the  total  energy 

JV  12 

functional  11  =  X]e=i  with  respect  to  the  displacements  Au^,  and  Au'^’’  respectively,  to  zero,  or 
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[  [{(Jij  +  Aaij)njSuf  dSfi  -  [  {ti  +  Ati)  ]  Suf  dT 
^  Jan.  Jr,„ 


=  0 


Tt. 

V  5uf  G  Vf  =  {vf  G  :  vf  =  0  on  r„e}  V  e  on  5fi 


(4.6) 


and 


J ^  [{(Tij  +  Aaij)nJ  -  cj)\]  5uf  dTor  =  0 

j ^  [{aij  +  Aaij)nJ  +  dTcr  =  0 

V  G  vr,  V  e  on  T,r  (4.7) 

cr  cr  Ir  Ir  12 

where  cj)  =  t’l°’^d{uf  -  u^)  is  the  cohesive  energy  function  and 

U'r^—U'r^  * 

4.2.1  Cohesive  zone  models  for  crack  propagation 

Cohesive  zone  models,  introduced  in  [5,  25]  and  developed  in  [63,  64,  65,  97,  30,  33,  39,  68],  are  effective  in 
depicting  material  failure  as  a  separation  process  across  an  extended  crack  tip  or  fracture  process  zone.  In 
these  models,  the  tractions  across  the  crack  reach  a  maximum,  subsequently  decrease  and  eventually  vanish 
with  increasing  separation  across  the  crack.  The  cohesive  model  used  in  this  chapter  is  a  three  parameter 
rate  independent  linear  cohesive  model,  proposed  in  [39,  68].  This  is  an  extrinsic  (two  stage)  model  which 
has  an  infinite  stiffness  or  slope  in  the  rising  portion  of  the  traction-separation  law  up  to  a  peak  traction 
value.  This  is  followed  by  linear  descending  segment  till  a  zero  traction  value  is  reached.  The  model  assumes 
a  free  cohesive  energy  potential  ^  such  that  the  traction  across  the  cohesive  surface  is  expressed  as: 


t 


coh 


(4.8) 
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Here  6„  and  dt  correspond  to  the  normal  and  tangential  components  of  the  opening  displacements  over  the 
cohesive  surface  in  the  n  and  t  directions  respectively.  An  effective  opening  displacement  is  defined  as 


s=y/sf^ 


(4.9) 


where  ^  is  a  coupling  coefficient  to  allow  assignment  of  different  weights  to  normal  and  tangential  opening 
displacements.  Consequently  the  cohesive  surface  traction  reduces  to 


^coh  ^  +  <5„n),  where  t=^  =  \J 


(4.10) 


where  and  are  the  normal  and  tangential  components  of  surface  tractions.  The  effective  cohesive 
force  t  in  this  model  for  increasing  d  takes  the  form 


t=< 


Vd<de 


Vd>  A 


(4.11) 


6e  corresponds  to  the  separation  at  which  t  goes  to  zero  and  cimax  is  the  peak  value  of  t.  The  effective  normal 
traction-separation  response  of  this  model  is  depicted  in  figure  (4.2).  In  the  softening  region,  Unloading  from 
any  point  on  the  traction-separation  curve,  proceeds  along  a  linear  path  from  the  current  position  to  the 
origin  as  shown  by  the  line  BO  in  figure  4.2.  The  corresponding  t  —  6  relation  is 


t  = 


A  -  Sr, 


Vd  <  Smax  <  d, 


(4.12) 


Reloading  follows  the  path  OBC  with  a  reduced  stiffness  in  comparison  with  the  original  stiffness.  Traction 
vanishes  for  d  >  dg. 

For  negative  normal  displacement  (compression),  stiff  penalty  springs  with  high  stiffness  are  introduced  be¬ 
tween  the  node-pairs  on  the  crack  face.  To  define  the  tangent  stiffness  matrix,  it  is  necessary  to  distinguish 
between  crack  initiation  (d  =  0)  and  crack  propagation  from  an  initialized  state  (d  >  0).  In  the  former. 
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=  t,  and  =  0  are  assumed,  which  implies  that  the  initiation  is  in  pure  mode  I.  The  cohesive  pa¬ 
rameters  in  this  study  are  calibrated  from  experiments  done  for  epoxy-steel  composites  as  discussed  in  [53,  35]. 

Recent  experimental-computational  studies  on  composites,  conducted  in  [93]  show  that  the  three  or  four 
parameter  cohesive  models  are  more  suitable  for  modeling  interfacial  debonding  in  comparison  with  the  two 
parameter  models  based  on  Ferrante’s  law  [63,  64,  65].  Similar  conclusions  have  also  been  drawn  in  the  work 
by  Ghosh  et.  al.  [35,  53],  where  bilinear  cohesive  models  were  chosen  to  study  interfacial  debonding  in  fiber 
reinforced  composites. 

4.2.2  General  element  assumptions  and  weak  form 

In  the  absence  of  body  forces,  two  dimensional  stress  fields  satisfying  equilibrium  relations  can  be  generated 
from  the  Airy’s  stress  function  $(x,y).  In  the  incremental  formulation,  stress  increments  are  obtained  from 
derivatives  of  the  stress  functions  A$(x,  y)  as: 


j 

/  2  \ 

1 

^(Jyy 

= 

dx^ 

=  [P(x,y)]{A/3} 

(4.13) 

1 

^O'xy  j 

j 

y  dxdy  j 

where  {A/3}  is  the  column  of  unknown  stress  increment  coefficients,  associated  with  the  stress  interpolation 
matrix  [P(x,y)].  Convergence  properties  and  efficiency  of  VCFEM  depend  on  the  choice  of  $.  These 
functions  should  adequately  account  for  the  geometry  and  location  of  the  heterogeneity  in  the  element. 
Polynomial  functions  alone  do  not  contribute  to  this  requirement  and  hence  lead  to  poor  convergence  [61,  80]. 
Consequently,  stress  functions  in  X- VCFEM  are  constructed  from  different  expansion  functions  that  have 
complementary  effects  on  the  solution  convergence  for  the  propagating  crack.  Compatible  displacement  fields 
satisfying  inter-element  continuity  on  the  element  boundary  d^lf  and  intra-element  continuity  on  the  crack 
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face  Ter  are  generated  by  interpolation  of  nodal  displacements,  [35,  61,  80]  as: 


{Aw"}  =  [Le]{Ag"}  on  Sfle 
{Ai"’-}  =  [W]{aJ"’'}  on  Ter, 

{A^’-}  =  [WjlAg"’-}  on  Tlr  (4.14) 

1  2 

The  interpolation  matrices  [Le],  [Lcr],  [Lcr]  for  the  nodal  displacements  on  the  respective  boundaries  are 
constructed  using  standard  linear  or  hierarchical  shape  functions. 

Remark  It  is  desirable  that  the  displacement  interpolations  on  the  crack  surface  in  equation  (4.14)  have  ad¬ 
equate  resolution,  consistent  with  the  high  resolution  in  the  stress  fields  near  the  crack  tip.  To  accommodate 
this,  hierarchical  shape  functions  are  added  to  standard  linear  shape  functions  to  describe  displacements  on 
the  crack  surface  as: 


4 

=  ^  7Vi(.)*qr  (4.15) 

i=l 

where  Ni  =  |(1  —  s),  N2  =  |(1  -I-  s),  N3  =  |(s^  —  1),  and  N4  =  |(s®  —  s).  The  first  two  are  the  standard 
linear  shape  functions,  while  the  last  two  are  the  hierarchical  shape  functions  in  natural  coordinates  s.  The 
degrees  of  freedom  corresponding  to  higher  order  shape  functions  (i.e.  to  quadratic,  cubic,  etc.)  cannot  be 
interpreted  as  nodal  values  of  displacement.  Instead,  they  are  values  of  some  higher  order  derivatives  of  the 
solution  at  the  midpoints  (or  linear  combination  of  these  derivatives). 

Substituting  the  interpolations  of  stress  and  displacement  fields  from  equations  (4.13)  and  (4.14)  into 
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equation  (4.2)  results  in  the  matrix  form  of  the  element  complimentary  energy 


where 


He  = 


+ 


+  A;9}  +  W  +  [G]'^{q'=  +  Aq-^} 

+  A;a}^[G=’-]{q-  +  Aq-}  -  +  Aq-^} 

2  2  2 

{^a  +  A;a}^[G=’’]{q“’’ +  Aq^"} 


indVcr 


[H]  = 

[  [P]^[S][P]dfl  [G1  = 

[  [P]"’[n'=][Le]d5fl 

'afie 

[G^n  = 

L  [Pf[n^^][L\r]dT„  [G“1  = 

/  [P]^[n-][Ll]dF, 

{t}  = 

(  {t  -I-  At}  \Le\dVtm 

(4.16) 


(4.17) 


Construction  of  appropriate  stress  functions  with  optimally  high  resolution  is  necessary  for  accurately  de¬ 
picting  high  stress  gradients  near  the  crack  tip. 


4.2.3  Stability  conditions 

Following  the  stability  conditions  derived  for  displacement-based  and  stress-based  finite  element  approxima¬ 
tions  in  [4,  12,  104],  the  stability  conditions  of  the  stress-displacement  field  variational  problem  in  X-VCFEM 
depend  on  the  following  conditions. 

•  The  matrix  [H]  should  be  positive  definite.  From  the  definition  of  [H]  in  equation  (4.17),  the  necessary 
condition  for  it  to  be  positive  definite  is  that  the  compliance  tensor  [S]  be  positive  definite,  which  is 
true  for  elastic  problems. 

•  A  second  condition  is  that  the  finite-dimensional  stress  subspaces  T  be  spanned  uniquely  by  the  basis 
functions  [P].  This  is  satisfied  by  assuring  linear  independence  of  the  columns  of  basis  functions  [P], 
which  also  guarantees  the  invertibility  of  [H]. 
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•  Additional  stability  conditions  should  be  satisfied  to  guarantee  non-zero  stress  parameters  /3  for  all 
non-rigid  body  displacement  fields  on  the  element  boundary  uf  or  on  the  crack  face  uj;’'.  This  is 
accomplished  by  careful  choice  of  the  dimensions  of  the  stress  and  displacement  subspaces,  i.e.  np  > 
rig  +nq''*2  —  3,  where  np  is  the  number  of  /3  parameters,  and  nf  and  are  the  number  of  displacement 
degrees  of  freedom  on  the  element  boundary  and  crack  face  respectively. 

4.3  Creation  of  Enriched  Stress  Functions  in  X-VCFEM 

VCFEM  formulations  for  micromechanical  analysis  of  heterogeneous  materials  have  incorporated  polynomial 
and  reciprocal  stress  functions  based  on  analytical  micromechanics  results  in  [35,  61,  80,  36].  In  the  present 
work,  the  heterogeneity  is  in  the  form  of  an  evolving  cohesive  crack.  Two  conditions  need  to  be  considered 
in  the  choice  of  stress  functions.  The  first  is  that  it  should  adequately  represent  crack  tip  high  stress 
concentration  as  required  by  the  cohesive  zone  models.  Polynomial  functions  alone  are  unable  to  satisfy 
this  requirement  and  hence  suffers  from  poor  convergence.  The  second  condition  is  that  the  stress  function 
should  account  for  stress  jump  across  the  crack  surface.  The  stress  functions  in  X-VCFEM  incorporate  three 
different  components,  namely:  (a)  a  purely  polynomial  function  to  yield  the  far  field  stress  distributions 
away  from  the  crack  tip,  (b)  a  branch  function  ^branch  jg  constructed  from  level  set  functions,  and  (c) 
a  multi-resolution  wavelet  function  to  account  for  the  moving  crack  tip  stress  concentration.  Thus, 

^  _  ^poly  _|_  ^branch  _|_  ^wvlt 


4.3.1  Pure  Polynomial  Forms  of  Stress  Function: 

The  pure  polynomial  component  of  the  stress  function  is  written  in  terms  of  scaled  local  coordinates 
)  V  =  with  origin  at  the  element  centroid  {xc,  Vc),  as: 

Pn,qn 

^po‘y(lfj)=  ^  (4.18) 

The  scaling  parameter  in  the  coordinate  representation  is  i  =  y^max(x  —  x^)  x  max(y  —  i/c) 

y{x,y)  G  dVtg.  The  use  of  the  scaled  local  coordinates  as  opposed  to  global  coordinates  {x,y)  in  the 
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construction  of  stress  functions,  prevents  ill  conditioning  of  the  [H]  matrix  due  to  the  high  exponents  of 
{x,y)  in  ^P°‘y.  As  discussed  in  [88],  invariance  of  stresses  with  respect  to  coordinate  transformations  can 
be  ensured  by  a  complete  polynomial  representation  of  ,  while  stability  of  the  algorithm  requires  linear 
independence  of  the  columns  of  stresses  derived  from  ^P°‘y. 

4.3.2  Branch  Stress  Functions  Using  Level  Set  Methods 

The  branch  function  ^branch  facilitates  jumps  in  stresses  across  the  crack  surfaces.  These  functions  should 
not  affect  the  solutions  in  the  continuous  region  beyond  the  crack.  This  construction  requires  a  functional 
representation  of  the  surface  or  line  of  discontinuity.  Level  set  methods,  introduced  by  Sethian  [1,  85]  for 
following  the  evolution  of  interfaces,  is  ideal  for  representing  arbitrary  contours.  The  method  has  been  used 
by  Belytschko  and  coworkers  in  [10]  for  the  construction  of  branch  functions  associated  with  the  partition  of 
unity  in  a  displacement  based  FEM  formulation.  The  standard  level  set  methods  invoke  continuous  evolution 
of  the  entire  surface  of  discontinuity.  However  for  problems  involving  cracks,  the  only  evolution  occurs  at  the 
crack  tip  and  the  crack  surface  needs  to  be  frozen  behind  tip.  A  vector  level  set  method  has  been  developed  in 
[100,  99]  to  freeze  the  crack  surface  in  accordance  with  geometric  updating.  This  method  is  used  in  this  work. 

An  approximation  to  the  crack  surface  Ter  in  figure  4.1  is  constructed  to  describe  the  discontinuous 
stress  fields  across  crack  paths.  As  shown  in  figure  4.3(a),  the  discontinuous  surface  is  expressed  by  a  signed 
distance  function  /(x)  defined  as 


/(x)  =  min  II  X  —  X  II  sign(n''“  •  (x  —  x))  (4.19) 

xGF 

where  x  is  a  point  on  the  surface  of  discontinuity  and  n+  is  a  unit  normal  pointing  in  the  direction  of  the 
region  of  positive  distance  function. 

Consequently,  x  is  the  closest  point  projection  of  any  point  x  on  Ter-  In  order  to  describe  the  crack 
path  accurately,  the  signed  function  /(x)  is  evaluated  at  every  integration  point  in  the  Voronoi  cell  element 
directly.  The  process  of  constructing  branch  functions  involves  steps  that  are  described  below. 
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•  Radial  distance  functions  to  the  two  crack  tips  ri  (x)  and  r2  (x)  and  the  corresponding  angular  positions 
01  (x)  and  02  (x)  are  depicted  in  figure  4.3(a).  These  functions  are  expressed  in  terms  of  coordinates  of 
local  systems  with  origins  at  the  crack  tips.  For  the  local  system  at  crack  tip  1,  the  coordinates 
of  X  are  (^i,7?i).  In  accordance  with  the  definition  of  the  signed  distance  function,  the  radial  distance 
and  angle  functions  are  expressed  as 


ri(x)  =  \J^l  +‘01  and  0i{x) 


TT  -  ^1  <  0,  /  >  0 

-  TT  6  <  0,  /  <  0 

6  >  0 


(4.20) 


Similarly,  the  radial  distance  and  angle  functions  for  the  coordinate  system  at  crack  tip  2  are  defined 
as: 


r2 


(x)  =  \J^l+ril 


and  02{x) 


TT  -  6  <  0,  /  >  0 

-  TT  6  <  0,  /  <  0 

sin-^  ^  6  >  0 


(4.21) 


•  The  branched  stress  function  is  constructed  in  terms  of  the  functions  /(x),  0i,  ri,  02,  and  r2,  as: 


^branch 


8=0, t=0 


(4.22) 


The  terms  rf  and  in  ^branch  necessary  for  avoiding  crack  tip  singularity  in  the  stresses  due  to 
this  function  and  for  improving  the  accuracy.  Along  the  tangential  extension  to  the  crack  path  at  the 
tip  1,  is  zero  since  sin^  =  0.  Hence  does  not  contribute  to  the  stresses  ahead  of  the 

crack  tip  1.  In  an  analogous  manner,  ^branch  ^ero  along  the  extension  to  the  crack  path  at 

the  tip  2,  since  cos^  =  0.  Therefore  does  not  contribute  to  the  stresses  in  this  region  also. 

However,  along  the  crack  surface  between  the  two  crack  tips,  sin^  =  ±1  on  both  sides  of  the  crack,  and 
cos^  =  1.  This  renders  in  equation  (4.22)  discontinuous  across  the  crack  path.  In  ^branch^ 

is  used  to  create  the  discontinuity  across  the  crack  surface,  while  02  eliminates  the  discontinuity  ahead 
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of  crack  tip  2.  In  some  special  instances  with  only  one  crack  tip  such  as  a  panel  with  an  edge  crack, 
equation  (4.22)  may  be  simplified  by  removing  r2  and  62  dependence  to  yield 


^branch  ^  ^ 

6,t  ^ 


(4.23) 


A  coordinate  transformation  is  required  to  obtain  stress  components  in  the  global  coordinate  system  from 
^branch  qjj  local  Coordinate  system. 


branch 


^  XX  ^  xy 
Cbxy  O'yy 


[Q6]^ 


Q2^branch 

Q2^branch 

d^idrii 


Q2^branch 

d^idrii 


Q2^branch 


Q^^branch 
q2  ^branch 

d^2  driz 


a_^ 


d^2  dr)2 


,2^branch 


~w 


[Q6] 


(4.24) 


where  [Qj]  is  the  transformation  matrix  from  (^i,7?i)  and  {^2,V2)  systems  to  {x,y),  and  is  expressed  as 


[Q6] 


dr}!  dr}! 

dy  dx 

dji  dji 

dy  dx 

dr]2  dr]2 

dy  dx 

de2  di2 

dy  dx 


(4.25) 


The  branch  function  is  evaluated  at  every  integration  point  in  the  element.  A  typical  function 
for  s  =  0  and  t  =  0  is  plotted  in  figure  4.3(b).  The  plot  shows  that  the  function  is  continuous  everywhere  in 
the  domain  except  across  the  crack  surface.  The  example  of  a  double  cantilever  beam  under  a  sliding  load, 
as  shown  in  figure  4.7,  explains  the  effect  of  level-set  method  based  branch  functions.  In  figure  4.7(a),  the 
dimension  is  a  =  1.5m.  Figure  4.7(b)  shows  the  stress  Uxx  plots  as  a  function  of  y  at  x  =  —0.3m.  The  stress 
functions  are  constructed  with  and  without  branch  functions  in  this  example.  Uxx  changes  its  sign  with  a 
jump  in  its  magnitude  on  different  sides  of  the  crack  and  the  jump  at  y  =  0  is  predicted  well.  However,  the 
transition  is  gradual  from  negative  to  positive  values  for  the  curve  without  branch  functions.  Although  the 
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transition  takes  place  in  a  short  interval,  the  method  is  not  able  to  catch  the  discontinuity  without  branch 

1  2 

functions.  This  also  results  in  the  matrices  [G^^]  and  [G^^]  in  equation  (4.16),  on  different  sides  of  the  crack 
to  be  linearly  dependent  on  each  other  (one  is  the  negative  of  the  other). 

4.3.3  Multi-resolution  Wavelet  Functions  for  Modeling  Cohesive  Cracks 

Wavelet  bases,  discussed  in  [16,  62],  are  L^iTZ)  and  generally  have  compact  support.  Only  the  local  coeffi¬ 
cients  in  wavelet  approximations  are  affected  by  abrupt  changes  in  the  solution,  such  as  for  shock  waves.  This 
localization  property  makes  the  wavelet  basis  a  desirable  tool  for  problems  with  a  high  solution  gradients, 
concentrations  or  even  singularity.  A  brief  introduction  to  wavelet  basis  functions  is  provided  next. 

Principles  of  wavelets  and  multi-resolution  analysis 

The  construction  of  wavelet  functions  starts  from  a  scaling  or  dilatation  function  and  a  set  of  related 
coefficients  {p{k)}kez  which  satisfy  the  two-scale  relation 

(p{x)  =  '^p{k)(p{2x  —  k)  (4.26) 

k 

The  scaling  function  has  a  compact  support  only  if  many  coefficients  p{k)  are  non-zero.  Translations  of  the 
scaling  function  ^(x  —  k)  form  an  unconditional  basis  of  a  subspace  Vq  C  L^{TZ).  Through  a  translation  of 
^  by  a  factor  of  2"  and  dilation  by  a  factor  of  k  ■  2^"  the  unconditional  basis  is  obtained  for  the  subspace 
Vn  C  L^{n)  as 

(x)  =2"/2^(2"x-fc)  (4.27) 

for  a  resolution  level  n.  The  scaling  function  (j)  is  defined  as  orthonormal  if  translations  at  the  same  level  of 
resolution  satisfies  the  condition 


/OO 

4>n,k{x)4>n,i{x)dx  =  5k,i  V  n,  fc,/ 

-OO 


G  .2 


(4.28) 
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Consequently,  the  best  approximation  of  a  function  f{x)  in  the  subspace  Vn  of  L^iTZ)  is  expressed  as  the 
orthogonal  projection  of  /  on  Vj,  as: 


/OO 

fix)4)„^kix)da 

-OO 


Approximation  of  f{x),  can  be  made  at  different  resolution  levels,  and  these  approximations  in  subspaces 
■  ■  ■  )  Ki-i!  Vn,  Ki+i!  ■  ■  ■ )  follow  the  relation 


{0}  =  y-oo  c  •  •  •  C  y-i  C  Vo  c  Vi  C  •  •  •  C  Voo  =  L^TZ),  where 
limn^<x>Vn  =  IJ  y^  is  dense  in  L‘^{TZ)  and  linin^-oo  n„  Vn  =  {0} 


In  the  multi-resolution  level  transition,  the  information  lost  in  the  transition  from  level  V„+i  to  level  V„  is 
characterized  by  an  orthogonal  complementary  subspace  Wn-  A  basis  for  the  subspace  Wn  can  be  obtained 
is  in  the  same  manner  as  for  scaling  function,  i.e.  by  dilating  and  translating  the  mother  wavelet  function 


V’(x)  =  E  q{k)'ip{2x  —  k) 


The  subspaces  spanned  by  the  wavelet  functions  have  the  following  essential  properties: 


(i)  V„+i  =  y,  0  W„  V,  i.e.  W„  is  the  orthogonal  complement  of  V„  toy,+i; 
{ii)  For  orthonormal  bases,  iy„i  is  orthogonal  to  iy„2; 

{in)  For  orthonormal  bases,  05^_oo  =  L‘^{TZ) 


An  approximation  of  the  function  f{x)  at  the  n  —  th  resolution  level  may  be  expressed  as  the  orthogonal 
projection  of  /  on  W„  as 


/OO 

f{x)'lpn,k{x)dx 

-OO 


48 


Due  to  the  orthonormality  and  multi-resolution  properties  of  wavelet  basis  functions,  higher  level  approximate 
solutions  can  be  generated  from  results  of  lower  level  solutions  (see  [16,  62])  by  selective  superposition  of 
complementary  solutions.  The  use  of  adaptive  enrichment  is  very  attractive  to  those  regions  where  a  pre¬ 
determined  ’error  or  residual’  tolerance  is  not  met  at  the  lower  level. 

Selection  of  the  wavelet  function 

Various  wavelet  functions  have  been  proposed  in  the  literature  for  numerical  solutions  of  ODEs  and  PDEs. 
These  functions  have  been  incorporated  in  the  method  of  weighted  residuals  like  the  Galerkin’s  method  and 
collocation  method  to  solve  problems  with  multi-level  features  in  [38,  47,  71].  Among  the  large  number  of 
wavelet  functions  proposed  are  the  Haar  function  [40],  the  Meyer’s  wavelets  [57],  the  Chui- Wang’s  B-spline 
wavelets  [17],  etc.  One  of  the  most  commonly  used  wavelet  functions  is  Daubechies’  compactly  supported 
orthonormal  wavelets  [21,  23,  38].  However,  they  are  constructed  through  recursive  algorithms  and  do  not 
have  an  explicit  analytic  expressions.  This  makes  it  is  difficult  to  obtain  their  first  and  second  derivatives, 
which  is  a  requirement  in  X-VCEEM  for  deriving  stresses  in  terms  of  stress  functions.  Also  the  orthonormality 
of  the  Daubechies  wavelet  cannot  be  transferred  to  the  orthonormality  for  stresses  by  differentiation,  and 
hence  they  are  not  considered  to  be  suitable  for  stress  functions  in  X-VCEEM.  Alternatively  a  family  of 
Gaussian  functions,  for  which  the  first  and  second  order  derivatives  are  popular  wavelets  bases  [11,  27,  52], 
is  implemented  in  the  representation  of  X-VCEEM  stress  functions  and  stresses.  The  expressions  for  the 
Gaussian  function  and  its  n  —  th  order  derivative  are: 

G'(x)  =  ^  (_l)n£L(e-(^)V2)  (4  34) 

The  dilation  and  translation  parameters  a  and  h  respectively  can  assume  arbitrary  values  and  can  be  changed 
in  a  continuous  fashion.  The  ability  of  wavelets  to  translate  diminishes  the  need  to  re-define  new  elements 
or  remesh  in  conventional  EEM  solution  of  problems  with  moving  boundaries.  By  changing  translation 
parameters,  the  multi-levels  of  wavelet  bases  can  be  made  to  closely  follow  a  moving  crack  tip.  Additionally 
the  dilation  parameter  with  compact  adjustable  window  support  can  be  used  to  provide  high  refinement 
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and  resolution.  Hence  it  is  a  convenient  way  of  moving  the  stress  concentrations  using  the  multi-resolution 
properties. 

Multiresolution  wavelet  based  stress  functions  for  crack  problems 

The  wavelet  based  stress  function  is  constructed  in  a  local  orthogonal  coordinate  system  (^,  rj) ,  centered  at 
the  crack  tip.  The  ^  direction  corresponds  to  the  local  tangent  to  the  crack  surface.  The  corresponding 
stress  function  $a,6,c,d  in  the  Gaussian  wavelet  basis  is  given  as: 


(4.35) 


where  a,b,c,d  are  parameters  that  can  take  arbitrary  continuous  values.  For  implementation  in  multi¬ 
resolution  analysis  involving  discrete  levels,  the  translation  and  dilation  parameters  should  be  expressed  as 
discrete  multiples  of  some  starting  values.  Consequently,  these  discrete  values  a™,  6„,  c*  and  di  are  expressed 
as: 

' 

am=ai- 
bn  —  Tl  '  b\  '  Clrn 

<  (4.36) 

Cfc  =  Cl  • 
di  —  I '  di  ' 

Here  {m,k)  correspond  to  the  levels  and  {n,l)  correspond  to  the  discrete  translation  of  the  bases  in  the 
{^,r})  directions  respectively.  The  parameters  (ai,ci)  are  the  initial  dilating  values  at  the  first  level  m  =  1, 
while  tra{<  1),  trc{<  1)  are  the  transfer  rates  from  one  level  to  the  next  higher  one.  The  parameters  bi,di 
represent  the  starting  values  of  a  step  translation  quantity  at  the  m  —  th  dilation  level.  The  narrow  (higher 
level)  wavelets  are  translated  by  small  steps,  whereas  the  wider  (lower  level)  wavelets  are  translated  by  large 
steps.  Parameters  tra  =  tfc  =  1  and  6i  =  di  =  0  imply  no  dilation  and  translation  respectively.  Parameters 
Co,  Cc,  and  do  are  counterparts  of  and  bo  in  rj  direction.  With  the  specific  relations  between  dilation 
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and  translation  parameters,  the  Gaussian  wavelet  enriched  stress  function  in  equation  (4.35)  becomes 


am 


A 


Tn,n,k,l 


(4.37) 


The  family  of  wavelet  enriched  stress  functions  in  equation  (4.37)  are  not  orthonormal,  but  they  construct 
a  linearly  independent  basis  [22].  This  leads  to  robustness  and  high  precision  in  the  reconstruction  of  any 
function  /  even  with  low  level  coefficients.  The  wavelet  enriched  stress  function  in  X-VCFEM  is  thus  written 
as 


m  u  I 

"’•n  5  2  i'“n 

m=l  ,n=—  ,k=l  ,1=0 


(4.38) 


The  corresponding  stresses  are: 


Q2^wvlt 
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[  didv 
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o  I  e  e  fc 
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(4.39) 


The  stress  components  in  the  global  coordinate  system  are  obtained  by  the  transformation  from  the  local 
coordinate  system  as 


p  n 
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(4.40) 


where  [Q.^,]  is  the  transformation  matrix  from  (^,7?)  to  {x,y): 
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(4.41) 


Figure  4.4  shows  the  support  region  for  the  wavelets  enriched  in  a  X-VC  element.  This  region  is 

positioned  symmetrically  in  the  vicinity  of  evolving  crack  tips.  The  crosses  (x)  corresponds  to  the  position 
of  each  wavelet  basis  function  d„  at  a  lower  level,  and  the  squares  (□)  correspond  to  additional  locations 
at  a  higher  level  in  the  multi-resolution  algorithm.  Only  the  points  at  the  top  half  are  shown  in  the  figure 
due  to  symmetry. 

The  method  of  implementation  of  the  multi-resolution  wavelet  enriched  stress  functions  in  X-VCFEM  is 
described  below. 

1.  For  the  starting  level  m  =  fc  =  1,  20  points  marked  by  crosses  (x)  in  figure  4.4  (a),  are  used  to  delineate 

the  wavelet  enriched  function  in  equation  (4.38).  This  corresponds  to  m  =  1,  n  =  5,  k  =  1 

and  I  =  4. 

2.  With  ensuing  higher  levels  in  the  multi-resolution  wavelet  functions  according  to  the  equation  (4.36), 
higher  level  wavelet  bases  are  added  to  the  stress  function  as  marked  by  squares  (□)  in  figure  4.4 
(b).  The  addition  is  done  adaptively  in  accordance  with  error  criteria  discussed  in  section  (4.3.4).  A 
refinement  in  the  starting  region  of  wavelet  enrichment  occurs  in  each  added  level,  i.e.  the  window  size 
of  additional  wavelet  basis  functions  is  smaller  than  ones  at  a  lower  level.  This  allows  a  zoom  in  to 
catch  higher  gradients  that  are  missed  at  the  coarser  scales. 

3.  The  process  of  successive  multi-level  refinement  can  continue  till  a  predetermined  error  tolerance  is 
reached. 

Remark  The  line  of  the  cohesive  crack  is  likely  to  intersect  the  region  of  support  of  the  wavelet  bases  func¬ 
tions.  It  is  important  for  the  numerical  algorithms  to  assure  that  wavelet  functions  based  on  one  side  of  the 
cohesive  crack  does  not  contribute  to  stresses  on  the  other  side.  The  influence  of  wavelet  stress  functions 
should  be  cut  off  across  this  line  of  discontinuity  by  establishing  a  truncated  effective  support  domain  for  the 
wavelet  function.  This  is  accommodated  by  ignoring  the  contribution  of  quadrature  points  in  the  numerical 
integration  on  the  other  side  of  the  crack  as  detailed  in  section  4.5.3. 
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In  summary,  the  stresses  in  an  element  are  computed  by  adding  contributions  from  equations  (4.20) 


(4.22)  and  (4.39),  to  yield 
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(4.42) 


4.3.4  Error  measure  for  adaptive  wavelet  enrichment 

The  Euler  equation  (4.3)  indicates  that  the  error  in  the  kinematic  equation,  which  is  satisfied  in  a  weak 
sense,  may  be  primarily  attributed  to  the  lack  of  adequate  resolution  in  the  equilibrated  stress  fields.  A 
strain  energy  based  element  error  measure,  derived  in  [80],  is  extended  to  the  present  problem.  Let  a  stress 
field  be  enriched  from  a  level  n  to  level  n  +  1  by  adding  the  wavelet-based  enrichment  stress  cr®"’’,  i.e. 


^level{n+l)  _  ^levelin)  _|_  ^enr 


(4.43) 


The  corresponding  percentage  change  in  the  strain  energy  {SE  =  cfijSijkicfki  dfi),  may  be  expressed  as 


A  SE  = 


5E(cr'®'®®'-"+i) 


xl00% 


(4.44) 


In  view  of  the  local  properties  of  wavelets  and  stress  concentration  at  crack  tips,  the  strain  energy  in 
equation  (4.44)  is  calculated  only  in  a  small  region  around  crack  tip  fienr  ■  Adding  levels  is  conditioned  upon 
the  requirement  that  ASE  is  less  than  a  preset  tolerance,  which  in  this  work  is  chosen  to  be  «  4%. 
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4.4  Solution  Method 


Crack  growth  in  multiply  cracked  materials  is  solved  using  an  incremental  approach,  where  a  set  of  elemental 
and  global  equations  are  solved  in  each  increment  for  stresses  and  displacements. 


1.  Local  equations  for  each  element  are  obtained  by  substituting  the  stress  interpolations  of  equation  (4.42) 
and  boundary /crack  face  displacement  interpolations  of  equation  (4.14)  in  the  element  energy  functional 
equation  (4.16)  and  setting  its  variation  with  respect  to  the  stress  coefficients  A;9  to  zero.  This  results  in 
the  weak  form  of  the  element  kinematic  relations 


[HUl3  +  A/3%  = 


[G«]  [G'n  -  [G" 


q®  +  Aq® 

q®’’  +  Aq" 

2  2 

qcr  ^  ^qcr 


(4.45) 


or  in  a  condensed  form 


[H]e{;a  +  A;a}e  =  [G]e{q  +  Aq}e 


(4.46) 


Since  equation  (4.46)  is  linear,  the  stress  coefficients  can  be  directly  expressed  in  terms  of  the  nodal  displace¬ 
ments,  provided  the  element  [H]e  matrix  is  invertible. 


2.  Subsequently,  the  weak  forms  of  the  global  traction  continuity  conditions  are  solved  by  setting  the 
variation  of  the  total  domain  energy  functional  with  respect  to  the  generalized  displacement  components  to 
zero.  This  results  in  the  weak  form  of  the  traction  reciprocity  conditions 
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54 


or  in  a  condensed  form: 
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The  forces  at  the  crack  surface  are  expressed  in  terms  of  the  cohesive  energy  as 
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Combining  equations  (4.46)  and  (4.48)  and  eliminating  the  stress  coefficients  {j3  +  Aj3}e,  results  in  the 
equation  for  solving  the  generalized  displacements 


JV 


JV 


E{[C^]nH]-MG]e}{q  + Aq}  =  E{Te.i}e 


(4.50) 


Equation  (4.50)  is  a  nonlinear  matrix  equation  system  due  to  the  cohesive  laws.  Consequently,  a  Newton- 
Raphson  iterative  solver  is  invoked  to  solve  for  the  increments  of  nodal  displacements.  The  linearized  form 
of  equation  (4.50)  for  the  j-th  iteration  is 


j (4-51) 

le=l  e=l  J 


which,  in  a  condensed  form  is 


[K^Ydq^  =  {RLJ  -  {RL}^' 


(4.52) 


A  numerical  problem  associated  with  modeling  cohesive  crack  growth  is  the  occurrence  of  snap-back  as  is 
shown  in  the  macroscopic  load-deformation  behavior  plot  of  figure  4.5. 

This  has  been  discussed  for  a  three  point  bending  solution  in  [59].  For  a  deformation  controlled  process 
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with  monotonically  changing  deformation,  the  solution  ignores  the  reverse  portion  of  the  displacement  BCD, 
occurring  with  snap-back.  The  Newton-Raphson  solver,  where  the  loading  process  is  monotonically  controlled 
by  incremental  deformation  or  load  conditions,  exhibits  a  discontinuous  drop  from  point  B  to  point  D.  It 
is  obvious,  that  this  solver  needs  to  be  augmented  with  the  capability  to  account  for  the  part  BCD,  i.e.  to 
decrease  both  load  and  deformation  with  the  growth  and  opening  of  the  crack.  The  arc-length  solver  has 
been  proposed  in  [19,  20,  84]  as  a  method  of  overcoming  this  shortcoming  by  introducing  an  unknown  loading 
parameter  (A  -I-  dX)  to  govern  the  load  increments.  Equation  (4.52)  is  modified  with  this  loading  parameter 
as 


[KS]^dq^'  =  (A^-  +  dA^){RLJ  -  {R-LF  (4.53) 

where  both  dX^  and  dq^  are  unknowns,  and  dX^  can  be  either  positive  or  negative.  The  additional  unknown 
dX^  requires  the  solution  of  a  constraint  equation,  written  in  terms  of  the  magnitude  of  the  deformation  of 
all  the  nodes  on  the  crack  surface  as 


^  ((A<)2  +  (A<)2)  =  AP  (4.54) 

iECrk 

where  Crk  represents  the  set  of  all  nodes  on  crack  surfaces.  A  summary  of  the  solution  process  is  explained 
in  the  flowchart  of  figure  4.6. 

4.5  Aspects  of  Numerical  Implementation 

4.5.1  Adaptive  criteria  for  cohesive  crack  growth 

A.  Direction  of  incremental  cohesive  crack  advance:  In  linear  elastic  fracture  mechanics,  it  is  common  to  use 
the  “maximum  hoop  stress  criterion”  to  determine  the  direction  of  crack  propagation  [7,  10].  Cracks  are 
assumed  to  propagate  in  a  direction  normal  to  the  maximum  hoop  stress  in  this  criterion.  Since  stresses  at 
crack  tip  are  singular  in  LEFM,  stress  intensity  factors  are  usually  used  to  determine  the  direction  of  crack 
propagation.  This  criterion  is  only  suitable  for  K-dominated  problems,  where  the  size  of  the  fracture  process 
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zone  is  small  compared  to  the  size  of  the  specimen.  A  different  criterion,  based  on  the  cohesive  energy  at 
the  crack  tip  is  used  in  X-VCFEM.  A  relation  between  the  cohesive  energy  (j)  for  complete  decohesion  and 
the  critical  energy  release  rate  Gc  has  been  established  in  [68]  from  the  definition  of  the  J— integral  as: 

35 

Gc  =  J=  t  — — dxi  =  /  td6  =  (j)  (4.55) 

Jo  dxi  Jo 

where  R  is  the  length  of  the  cohesive  zone.  Consequently,  the  crack  growth  direction  is  estimated  as  that, 
along  which  Gc  or  equivalently  the  cohesive  energy  (j)  is  maximized  for  a  given  crack  tip  state  of  stress.  The 
cohesive  energy  (j)A  at  the  crack  tip  A  along  any  direction  a  can  be  expressed  for  an  arbitrary  separation 
d(a)  as: 

a)  \  /  /■*(«)  Q§  \ 

t{a)d6j  ~  yj  y  (4.56) 

where  t{a)  =  is  the  magnitudes  of  the  effective  cohesive  force.  The  corresponding 

unit  normal  n  and  tangential  t  vectors  along  the  direction  a  are  expressed  as 

n  =  — sinai  +  cosaj  ,  t  =  cosai  +  sinaj  (4-57) 


(pAia)  = 
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The  normal  and  tangential  components  of  the  cohesive  traction  force  at  an  angle  a  may  then  be  deduced  as: 
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and  hence  the  effective  cohesive  traction  for  direction  a  is 


t{a)  = 

^((TssSin^a  —  CTa;ySin2a  +  GyyCos^aY  +  -(Ta;a;Sin2a  +  axyCos2a  +  -ayysm2ay 


The  incremental  direction  of  crack  propagation  is  assumed  as  that  which  maximizes  the  cohesive  energy  at 
A,  according  to  the  criterion: 


2^=0  and 

aa  aa^ 


A  combination  of  equations  (4.56),  (4.59)  and  (4.60),  yield 


</-^(a)  =  ^(aL.-i(«)^) 

^^m.ax 


d(j)A  _  ^  dt  _  ^ 

dot  ^max  dot 


—  =  [(CTa;a;Sin^a  —  axyS\.n2a  +  CTyyCOS^a)(CTa;a;Sin2a  —  2axyCOs2a  —  ayySin2a)  + 

/3^^(— ^CTa,3,sin2a  +  CTa,j^cos2a  +  ^CT3,a,sin2a)(— CTa,3,cos2a  —  2CTa,j^sin2a  +  CTj^j^cos2a)]/ 

1  J  0.5 

((CT3,a,sin^a  —  CTa,j^sin2a  +  CTj^j^cos^a)^  +  /3  ^(— -CT3,a,sin2a  +  CT3,j^cos2a  +  -CTj^j^sin2a)^) 


d‘^(l)A  Se 


-\{(JxyAB.2a  —  CTssSin^a  —  CTyyCos^a)(CTa;a:Sin2a  —  2CTa;yCos2a  —  CTyySin2a) 


+/3  ^(-CTa;a;Sin2a  —  CTa;yCos2a  —  -CTyySin2a)^(— CTa;a;Cos2a  —  2axyAri2a  +  CTyyCos2a)] 


Equation  (4.61)b  results  from  the  fact  that  t  cannot  be  equal  to  zero  for  decohesion  to  initiate  and  hence 
the  necessary  condition  evolves  from  its  derivative.  The  direction  of  crack  propagation  etc  is  obtained  as  the 


solution  of  equation  (4.61  )b  as 
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(4.62) 


The  optimal  angle  a^vcFEM  jg  chosen  as  the  one  that  satisfies  the  condition  in  equation  (4.61)c.  The 
corresponding  angle  given  by  the  maximum  hoop  stress  criterion  in  LEFM  is  expressed  in  terms  of  the  stress 
intensity  factors  Kj,  Ku  as: 
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(Jtc 
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{Ki/Kii  ±  ^/{Ki/Kuf  +  s) 


(4.63) 


where  the  sign  is  chosen  to  make  the  hoop  stress  positive.  The  first  of  equation  (4.62),  which  is  the  only 
choice  for  /?  =  1,  exactly  matches  the  angle  given  by  the  maximum  hoop  stress  criterion  (4.63).  For  the  pure 
sliding  problem  shown  in  figure  4.7,  etc  predicted  by  the  equation  (4.63)  is  70.5°,  while  that  by  X-VCFEM 
for  cohesive  stresses  is  68.2°. 


B.  Length  of  the  incremental  cohesive  crack  advance:  Upon  establishing  the  direction  of  incremental  cohe¬ 
sive  crack  growth  Uc,  the  length  of  cohesive  zone  advance  (A/)  should  be  estimated  in  the  crack  evolution 
scheme.  The  criterion  used  is  that  the  cohesive  energy  goes  zero  at  the  end  of  the  new  segment  as  shown  in 
figure  4.8(a).  To  achieve  this,  the  cohesive  energy  at  two  points  A  (present  crack  tip)  and  B  (close  to  A  in 
the  direction  of  crack  propagation)  are  evaluated  by  substituting  the  stresses  in  equation  (4.61)a.  The  tip  of 
the  cohesive  zone  is  obtained  from  the  linear  extrapolation  of  this  line  to  yield  zero  cohesive  energy.  From 
figure  4.8  (a),  the  increment  of  cohesive  crack  length  is  defined  as: 

A[=— (4.64) 

<PA  -  (pB 
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C.  Cracks  crossing  the  interelement  boundaries  and  merging  with  each  other:  Crack  advance  from  one  Voronoi 
cell  element  to  the  next  is  conducted  in  X-VCFEM  using  an  algorithm  depicted  in  figure  4.8(b).  A  contin¬ 
uous  tracking  method  is  implemented  to  monitor  if  a  cohesive  surface  has  reached  or  gone  past  an  element 
boundary.  In  this  method,  the  intersection  of  the  crack  surface  and  an  element  boundary  is  obtained  by 
solving  the  equation  system 

^  _  y  Vi  ^  _  y  yn 

^i+1  yi+1  yi  ^n+1  2/n+l  yn 

where  {xi,yi)  represents  the  tip  of  the  cohesive  crack  line  for  the  ith  increment,  and  {x„,y„)  is  the  position 
of  the  nth  node  on  the  element  boundary.  If  the  intersection  point  is  outside  of  the  cohesive  line  or  the 
element  boundary,  no  intersection  is  assumed.  Once  a  cohesive  crack  has  reached  its  intersection  with  the 
boundary,  a  new  node  pair  (m  n2)  is  introduced  on  the  element  boundary  at  this  point.  The  node  pair 
belongs  to  the  intersection  of  the  element  boundary  and  the  cohesive  crack,  i.e.  nin2  G  H  Fc^.  The 
crack  is  subsequently  advanced  to  the  next  element  following  the  usual  procedure  outlined  before. 

Another  condition  that  is  considered  in  this  work  is  the  merger  of  multiple  cracks  as  shown  in  figure  4.8(c) 
for  two  cohesive  cracks.  The  algorithm  for  crack  merging  is  an  extension  of  the  intersection  algorithm, 
discussed  above.  At  the  end  of  every  increment,  all  the  cracks  that  have  propagated  in  that  increment  are 
recorded.  Subsequently,  the  intersection  of  the  last  incremental  segment  of  the  cohesive  crack  with  those  of  all 
neighboring  cracks,  that  belong  to  either  the  same  element  or  neighboring  elements,  is  checked  using  equation 
(4.65).  Once  the  intersection  of  two  crack  segments  is  ascertained,  a  three-node  junction  (mi, m2,  m3),  as 
shown  in  figure  4.8(c),  is  inserted  at  the  point  of  intersection.  The  contribution  of  the  junctions  nodes  e.g. 
(mi,  m2)  to  the  load  vector  in  the  assembled  matrix  equation,  requires  special  treatment.  For  each  of  these 
nodes,  contributions  of  integrals  from  adjoining  crack  segments  belonging  to  two  different  cohesive  cracks, 
are  summed. 


(4.65) 
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4.5.2  Evaluation  of  stress  intensity  factors 


The  stress  intensity  factors  and  J—  integral  are  evaluated  in  the  post-processing  phase  of  the  computations. 
From  linear  fracture  mechanics,  the  relation  between  J—  integral,  stresses  and  stress  intensity  factors  are 
given  as 


-fA 


_Kj 


^ii 


E* 


(4.66) 


where  E*  =  E  (Young’s  modulus)  for  plane  stress,  E*  =  for  plane  strain,  and  p  is  the  Poisson’s 

ratio.  In  displacement  based  FEM  [59,  26],  the  contour  integral  is  converted  into  a  domain  integral  to 
improve  the  accuracy  of  the  stress  intensity  factors,  since  the  stresses  are  more  accurate  in  the  interior  of  an 
element.  However  in  X-VCFEM,  stresses  on  the  contour  and  the  interior  are  equally  accurate  due  to  stress 
interpolation  and  the  contour  integral  can  provide  similar  accuracy  as  the  domain  integral.  A  method  to 
extract  the  stress  intensity  factors  Kj  and  Ku  from  the  J-integral,  proposed  in  Yau  [105],  is  implemented 
in  X-VCEEM.  Displacement  fields  are  not  interpolated  in  the  interior  of  the  Voronoi  cell  element,  and  hence 
the  term  U2,i  in  equation  (4.66)  requires  a  special  evaluation  method. 

1.  Compute  eii,  €22,  and  ei2  at  a  series  of  points  {xi,yi,  i  =  in  a  small  shadowed  region  around 

the  integration  point  {xo,yo)  in  figure  4.8(d).  The  displacement  gradient  ui^i  is  calculated  from  en. 

2.  Eor  evaluating  W2,i  displacements  ui  and  U2  at  any  point  {xi,yi)  are  interpolated  using  polynomial 
functions. 


ui{xi,yi)  —  ao  +  aiXi  +  a2yi  +  asx^  +  ■  ■  ■  , 

U2{xi,yi)  =  ha  +  hiXi  +  h2yi  +  hzxl -  (4.67) 

where  Oq,  ai,  •  •  • ,  Om  and  bo,  h,  ■  ■  ■ ,  bu  are  unknown  coefficients.  To  constrain  the  rigid  body  motion, 
coefficients  are  evaluated  from  displacement  values  at  two  points  on  boundaries. 

3.  Displacement  gradient  expressions,  ui^i{xi,yi),  U2,2{xi,yi),  ui^2{xi,yi)  and  U2,i{xi,yi)  are  obtained  by 
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taking  derivatives  of  the  expressions  in  equation  (4.67).  Strain  expressions  in  terms  of  the  unknown 
coefficients  are  computed  from  these  derivatives.  At  each  point  {xi,yi),  the  strains  can  also  be  computed 
from  the  known  stresses  and  the  compliance  tensor,  i.e.  {e}  =  [S][P]  {j3}.  The  unknown  coefficients  oq, 
di,  ■  ■  ■ ,  dM  and  bo,  bi,  ■  ■  ■ ,  bu  in  equation  (4.67)  are  estimated  by  solving  a  least  square  minimization 
problem  for  the  strains.  Subsequently  the  displacement  gradient  U2,i  is  determined  at  the  integration 
point  {xo,yo)- 


4.5.3  Numerical  integration  schemes  for  [H]  and  [G]  matrices 

Integration  of  [H]  matrix: 

Numerical  integration  over  each  element  is  conducted  by  the  Gaussian  quadrature  method  to  form  the  matrix 
[H]  in  equation  (4.17).  In  this  method,  each  Voronoi  cell  element  is  recursively  subdivided  into  triangular 
subdomains,  on  which,  integration  points  are  generated  for  the  Gaussian  quadrature.  The  steps  involved  are 
discussed  below. 


1.  For  each  Voronoi  cell  element  shown  in  figure  4.9,  the  centroid  O  is  first  generated.  The  first  set  of 
triangular  subdomains  is  created  by  joining  each  of  the  vertices  of  the  cell  e.g.  (A,  B,  C,  D,  E, 
F)  with  the  centroid  O. 

2.  Each  triangle  is  further  subdivided  into  two  triangles  if: 


Area  of  triangle 
Area  of  Voronoi  cell  element 


>  TOLa 


(4.68) 


For  the  subdomain  triangle  BCO  shown  in  figure  4.9,  two  triangles  are  created  by  bisecting  the  longest 
edge  BC  at  O’  and  joining  it  with  the  opposite  vertex  O. These  new  smaller  triangles  are  again  checked 
against  the  tolerance  condition  and  further  dissection  is  executed  if  necessary.  Numerical  integration 
in  each  triangular  subdomain  is  done  using  13  Gauss  points. 

3.  For  the  region  containing  the  crack  tip  shown  in  figure  4.9,  a  smaller  value  of  TOLarea  is  chosen  in 
comparison  with  other  regions.  This  facilitates  a  higher  density  of  integration  points  in  regions  of  high 
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stress  gradients.  The  tolerance  in  an  element  is  consequently  adjusted  according  to  the  distance  of  the 
center  of  the  triangular  subdomain  from  the  crack  tip,  i.e. 

TOLarea  =  TOL^fa  +  ^  - r 

Jj 

where  i  is  a  scaling  parameter  defined  in  subsection  (4.3.1),  dtri  is  the  distance  of  the  crack  tip  from 
the  subdomain  and  ,  TOT™™  are  assumed  tolerances.  In  this  work  the  tolerances  are  chosen 

as  =  10%  and  TOL™™  =1%. 

4.  The  intersection  of  the  support  of  wavelet  functions  with  the  cohesive  crack  line  call  for  a  truncated 
support.  This  is  done  by  eliminating  the  contribution  of  quadrature  points  that  lie  on  the  other  side 
of  crack  face  from  the  wavelet  center.  A  visibility  criterion  introduced  in  [8]  provides  an  easy  way  to 
accommodate  this  discontinuity  in  the  construction  of  truncated  support.  In  this  method,  the  cracks 
are  considered  to  be  opaque  when  generating  valid  numerical  integration  regions.  A  ray  is  emitted 
from  the  center  W  of  a  wavelet  basis  function  in  an  arbitrary  direction  as  shown  in  figure  4.4.  If  it 
encounters  an  internal  crack,  the  ray  is  terminated.  All  quadrature  points  lying  in  the  dark  shadow 
region  on  the  other  side  of  the  crack  CC  are  suppressed  during  numerical  integration  of  this  wavelet 
basis. 


Integration  of  the  [G]  matrices: 

1  2 

In  equation  (4.17),  the  matrices  and  [G'^’’]  are  numerically  integrated  over  the  crack  surfaces  and  the 
matrix  [G1  over  the  element  boundary.  All  numerical  integrations  on  the  element  boundary  and  crack 
surfaces  are  executed  using  the  Gaussian  quadrature  method.  The  number  of  integration  points  Nint  on 
each  boundary /crack-face  segment  depends  on  the  distance  dsMe  between  its  center  and  the  crack  tip,  and 
is  chosen  from  the  condition 


9 


^int  —  \ 


dside  ^  O.IT 


16  dside  <  O.li 


(4.70) 
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where  L  is  the  scaling  parameter. 


4.5.4  Invertibility  of  the  [H]  matrix 

A  nonsingular  or  invertible  [H]  matrix  necessitates  the  linear  independence  of  the  columns  of  the  [P]  matrix. 
For  pure  polynomial  expansions  of  the  stress  functions,  this  condition  is  naturally  attained.  However  when 
adding  the  other  terms,  some  of  the  terms  in  the  branch  and  wavelet  functions  may  have  linear  dependence 
on  the  polynomial  terms.  In  X-VCFEM,  the  rank  of  the  [P]  matrix  is  first  determined  from  the  diagonal 
matrix  resulting  from  a  Cholesky  factorization  of  the  square  matrix 

[H*]  =  [  [P]^[P]dO  (4.71) 

Nearly  dependent  columns  of  [P]  will  result  in  very  small  pivots  during  Cholesky  factorization.  The  cor¬ 
responding  branch  and  wavelet  function  terms  are  dropped  from  the  stress  function  to  prevent  numerical 
inaccuracies  in  inverting  [H]. 


4.5.5  Elimination  of  element  rigid  body  modes 

X-VCFEM  uses  a  stress-based  formulation  with  independent  representation  of  displacement  fields  on  the 

element  and  crack  boundaries.  In  general,  the  nodes  of  the  crack  face  are  not  topologically  connected  to 

the  element  boundary  nodes.  However  it  is  important  that  all  nodes  in  the  element  possess  the  same  rigid 

body  modes.  The  rigid  body  modes  of  the  element  boundary  displacements  {q®}  are  directly  constrained 

in  the  solution  process  through  prescribed  displacement  boundary  conditions.  However,  it  is  necessary  to 

1  2 

connect  these  with  rigid-body  modes  for  the  crack  face  displacement  fields  {q®’’}  and  {q®'’}.  Singular  value 
decomposition  or  SVD  has  been  discussed  in  [80]  as  an  effective  method  for  identifying  and  constraining 
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rigid  body  modes  at  interfaces  inside  the  Voronoi  cell  elements.  The  matrix  product  may  be  expressed  as 

/  \ 

1 

=  [U][A][V]|  \  [  =  [U][A] 

qcr- 

<  J 

1 

=  [G"]  - 

r  1 

[U]  and  [V]  are  orthonormal  matrices  obtained  by  SVD  of  —  [G'^’’]  .  [A]  is  a  rectangular  matrix 

with  nonnegative  values  on  the  diagonal.  The  zero  or  singular  (very  small  values  in  numerical  computations) 
values  in  [A]  corresponds  to  either  trivial  solutions  or  rigid  body  modes  of  the  displacement  solution.  For 
accurate  displacements,  elements  in  corresponding  to  small  or  zero  eigen- values  in  [A]  are  eliminated. 


4.6  Numerical  Examples 

The  numerical  examples  solved,  are  divided  into  four  categories.  In  the  first  set  of  examples,  the  convergence 
of  X-VCFEM  enriched  by  multi-resolution  wavelet  functions  is  demonstrated  for  static  cracks  by  comparison 
with  theoretical  predictions  and  results  available  in  the  literature.  The  second  set  of  examples  show  the 
effectiveness  of  X-VCFEM  in  modeling  the  propagation  of  multiple  cohesive  cracks.  The  third  set  of  examples 
is  intended  to  investigate  the  effect  of  cohesive  parameters  on  crack  growth.  The  final  set  of  examples  looks 
into  the  growth  of  multiple  pre-existing  cracks  to  comprehend  the  effect  of  morphology,  e.g.  distribution, 
orientation  etc.. 

4.6.1  Convergence  tests  for  X-VCFEM  for  static  cracks 

Effects  of  translation  and  dilation  parameters 

Eigure  4.10(a)  shows  a  center  cracked  plate  of  width  2w=4cm  and  length  b=12cm  with  a  crack  length 
of  2a=1.6cm.  The  plate  is  loaded  in  simple  tension  with  a  constant  remote  load  of  cr  =  5  MPa.  The 
material  parameters  are:  Young’s  modulus  E  =  1  MPa  and  Poisson  ratio  p  =  0.3.  Due  to  problem 
symmetry,  only  the  right  half  of  the  plate  is  modeled  with  one  X-VCEEM  element,  as  shown  in  figure 
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4.10(b).  Symmetry  conditions  are  imposed  on  the  left  edge.  The  crack  face  is  modeled  using  10  node-pairs 
and  the  element  boundary  consists  of  22  segments.  The  stress  function  used  in  this  example  consists  of 
the  three  parts  discussed  in  section  4.3.  For  the  polynomial  function,  the  order  of  interpolation  in  equation 
(4.18)  corresponds  to  Pn  =  13  and  q„  =  13  for  a  total  of  102  terms.  For  the  branch  function  in  equation 
(4.22)  consists  of  only  1  term  with  =  0  and  =  0  .  The  wavelet  functions  are  changed  from  a  lower  level 
to  a  higher  level  using  the  adaptation  criterion  discussed  in  section  4.3.4.  Similar  parameters  are  assumed 
for  the  and  rji  directions,  i.e.,  ai  =  ci,  tra  =  tVc  and  hi  =  di  in  equation  (4.36).  The  starting  values  of 
the  parameters  for  the  lower  level  (m  =  fc  =  1)  are:  — 2  <  n  <  3,  0  <  /  <  1,  and  oi  =  ci  =  0.1.  The 
result  of  X-VCFEM  for  this  problem  is  plotted  in  terms  of  the  stress  Uyy  along  the  crack  face  {y  =  0) 
as  a  function  of  the  distance  from  the  center  of  the  crack  in  figure  4.11.  The  figure  4.11(a)  corresponds  to 
the  stresses  by  varying  the  translation  parameter  6i,  while  figure  4.11(b)  is  for  the  variation  of  the  dilation 
parameter  ai.  From  figure  4.11(a)  it  is  evident  that  a  smaller  hi  make  the  stress  concentration  at  the  crack 
tip  higher.  However,  very  small  hi  <  0.001  (no  translation)  leads  to  linear  dependence  of  the  columns  of 
the  [P]  matrix  generated  from  the  wavelet  basis  functions,  and  should  be  avoided.  Figure  4.11(b)  shows 
that  smaller  m  results  in  faster  convergence  to  higher  crack  tip  stress  concentration.  However,  very  small 
values  of  oi  can  also  lead  to  oscillatory  stresses.  On  the  other  hand,  large  ai  values  («  0.15)  shifts  the 
stress  peak.  The  optimal  selection  of  these  parameters  is  therefore  very  important.  This  is  obtained  through 
the  multi-resolution  construction  of  bases,  discuss  next.  The  multi-resolution  wavelet  bases  are  significantly 
more  effective  in  simulating  crack  problems.  Table  4.7  shows  the  effect  of  the  dilation  transfer  rate  tra  =  tfc 
on  the  stress  intensity  factors  for  variation  in  the  translation  parameters  hi  =  di.  Other  parameters  in  the 
simulation  are  oi  =  ci  =  0.1,  m„  =  k„  =  2>,  n„  =  6,  and  In  =  2.  The  values  tra  =trc  =  l  imply  no  dilation. 
As  tra  approaches  1,  the  different  levels  functions  become  more  and  more  dependent  on  each  other.  From 
the  table,  the  minimum  error  is  achieved  for  hi  =  0.1  and  tra  =  trc  =  0.5  or  0.6. 

Convergence  with  multi-resolution  wavelet  bases 

The  example  in  section  4.6.1  is  considered  again  for  studying  the  solution  convergence  behavior  with  multi¬ 
resolution  wavelet  functions.  The  four  sets  of  parameters  represent  four  instances  of  multi-resolution  stress 
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function  enrichment.  The  first  case  consists  of  only  polynomial  and  branch  functions  for  the  stress  inter¬ 
polation,  for  which  the  details  are  provided  in  section  4.6.1.  Cases  2,  3,  4,  and  5,  introduce  different  levels 
of  the  wavelet  basis  functions.  The  wavelet  parameters  common  to  these  four  cases  are:  n„  =  6,  /„  =  2, 
ui  =  Cl  =  6i  =  di  =  0.1,  and  tra  =  tVc  =  0.5.  The  parameters  corresponding  to  the  levels  of  the  multi¬ 
resolution  enrichment  (m„  =  fc„)  are  listed  in  table  4.7. 

The  mode  I  stress  intensity  factor  is  calculated  for  all  the  five  cases  and  is  normalized  with  respect  to  the 
analytical  prediction  K^-ef  by  linear  elastic  fracture  mechanics  (LEFM),  reported  in  [91].  The  second  row  of 
table  4.7  compares  this  value  for  the  different  cases.  Without  the  wavelet  bases  (case  1),  the  solution  is  16% 
higher  than  the  theoretical  value.  Cases  2-5  results  demonstrate  that  the  wavelet  basis  effectively  reduces 
the  error  with  increasing  resolution  level  (m„).  The  X-VCFEM  generated  stress  Oyy  at  y  =  0  is  plotted  in 
figure  4.12  for  cases  1-4.  Without  the  wavelet  enrichment,  the  stress  concentration  at  the  crack  tip  (x  =  0.8) 
is  completely  misrepresented.  The  stress  peaks  are  represented  with  increasing  accuracy  with  additional 
levels  of  multi-resolution  wavelet  functions.  The  strain  energy  error  in  equation  (4.44)  is  also  calculated  for 
the  cases  2-5  and  tabulated  in  table  4.7.  The  error  rapidly  decreases  with  increasing  wavelet  enrichment, 
confirming  the  fast  convergence  rate  of  the  multi- resolution  algorithm.  However,  the  stress  intensity  factor 
Ki  is  calculated  from  a  contour  that  is  away  from  the  crack  tip.  The  stresses  on  this  contour  are  much 
more  stabilized  and  additional  wavelet  bases  do  not  affect  these  stresses  considerably.  Hence,  the  error  in 
Kj  is  not  significantly  affected  by  their  addition.  From  the  above  convergence  tests,  the  optimal  parameters 
for  stress  function  representations  in  X-VCFEM  are  chosen  to  be  Pn  =  Qn  =  13,  Sn  =  tn  =  0,  n„  =  6, 
In  =  1,  m„  =  =  4,  ui  =  Cl  =  6i  =  di  =  0.1  and  tra  =  tfc  =  0.5.  These  are  retained  for  all  subsequent 

simulations.  X-VCFEM  simulations  of  the  cracked  plate  are  further  conducted  for  different  crack  lengths, 
to  study  the  effect  of  this  length  on  the  solution  convergence.  The  specific  dimensions  in  figure  4.10(a) 
are  2w=2  cm,  b=6  cm,  while  the  crack  length  2a  is  varied.  The  plate  is  loaded  under  remote  tension  of 
(T  =  40Pa.  X-VCFEM  solution  of  Ki  for  various  values  of  a/w  are  plotted  in  figure  4.13  and  compared 
with  the  theoretical  predictions  of  [91].  X-VCFEM  predictions  match  the  theoretical  results  extremely  well. 
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4.6.2  Efficiency  and  Accuracy  of  X-VCFEM 


Prior  to  studying  the  effect  of  cohesive  parameters  and  multi-crack  distributions,  the  accuracy  and  efficiency 
of  X-VCFE  model  are  validated  by  several  numerical  examples. 

Comparing  efficiency  with  ABAQUS  for  a  simple  crack  propagation  problem 

A  plate  with  a  pre-existing  edge  crack  under  remote  tension  load  is  solved  for  plane  strain  by  X-VCFEM 
and  ABAQUS  as  shown  in  figure  4.14(a).  The  material  Young’s  modulus  E  =  70,000MPa,  and  Poisson 
ratio  v  =  0.33.  A  bilinear  cohesive  zone  model  discussed  in  [53]  is  used  to  describe  the  crack  growth  and 
the  cohesive  model  parameters  are  amax  =  5  MPa,  dc  =  1  x  10^®  mm,  6e  =  5  x  10^®  mm  and  /3  =  0.707. 
The  entire  domain  is  represented  by  a  single  element  in  X-VCFEM,  consisting  of  142  nodes  for  displacement 
interpolation.  The  adaptive  enrichment  of  wavelet  bases  is  determined  by  the  strain  energy  error  in  equation 
(4.44).  As  shown  in  previous  section,  the  optimal  parameters  for  stress  function  representations  in  X-VCFEM 
are  chosen  to  be  Pn  =  Qn  =  13,  Sn  =  tn  =  0,  n„  =  8,  /„  =  1,  m„  =  =  4,  ui  =  ci  =  6i  =  di  =  0.1  and 

tTa  =  tVc  =  0.5,  which  means  that  stress  function  interpolations  consist  of  102  terms  of  polynomial  functions, 
1  term  in  the  branch  function,  and  128  terms  in  the  wavelet  function  representation.  These  are  retained  for 
all  subsequent  simulations.  It  is  assumed  that  the  crack  propagates  horizontally  due  to  symmetry,  and  hence 
the  modules  for  determining  incremental  crack  direction  in  section  4.5  is  switched  off  for  this  problem.  A 
special  UEL  subroutine  is  developed  in  ABAQUS  for  incorporating  the  cohesive  model  at  a  given  interfaces. 
A  total  of  12840  4-node  2D  element  and  77  cohesive  elements  are  used  in  ABAQUS.  Figure  4.14(b)  shows 
the  load  u-vertical  displacement  Uy  plot  at  point  A.  The  two  codes  yield  very  similar  results,  an  attestation 
of  X-VCFEM  accuracy.  However,  the  X-VCFEM  simulation  takes  only  1.6  minutes  on  a  single  CPU  in  the 
Pentium  4  cluster  with  2.4Ghz  Intel  P4  Xeon  processors,  as  opposed  to  13.9  minutes  by  the  ABAQUS  run  on 
the  same  machine.  Thus,  even  for  this  simple  example,  a  tenfold  advantage  in  computing  speed  is  achieved 
by  X-VCFEM.  It  is  expected  that  this  factor  will  increase  considerably  with  increasing  complexity,  such  as 
more  cracks. 
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A  classical  problem  on  dynamic  crack  propagation 


This  numerical  example  is  based  on  Kalthoff’s  well  known  experiment  on  dynamic  crack  propagation  in  a 
impact  loaded  prenotched  plate,  that  has  been  the  subject  of  many  studies  [50,  51,  76].  These  studies  suggest 
that  a  crack,  subjected  to  a  tension-compression  load  as  shown  in  figure  4.15(a),  propagates  at  an  angle  of 
approximately  60°  —  70°  with  respect  to  the  initial  notch  in  the  plate.  The  present  X-VCFEM  does  not 
explicitly  incorporate  inertia  terms,  and  hence  a  quasi-static  crack  propagation  problem  is  simulated  instead 
of  the  dynamic  test.  The  configuration  in  figure  4.15(a),  shows  that  the  experimental  projectile  motion  is 
replaced  by  the  traction  boundary  conditions  in  the  simulation  under  plane  strain  conditions.  A  small  initial 
crack  length  of  a=0.02m  is  chosen  to  mitigate  the  effect  of  the  constrained  right  hand  boundary  on  crack 
propagation.  Material  properties  for  this  problem  are:  Young’s  modulus  E  =  207  GPa  and  Poisson  ratio 
v  =  0.3  and  cohesive  zone  model  parameters  in  equation  (4.11)  are:  amax  =  0.1  MPa,  dg  =  1  x  10^®m,  and 
/3  =  0.  The  entire  domain  is  represented  by  a  single  element  in  X-VCFEM,  consisting  of  132  nodal  degrees 
of  freedom.  The  results  of  the  X-VCFEM  simulation  is  shown  in  figure  4.15.  From  figure  4.15(a)  the  initial 
crack  growth  angle  is  around  70°,  which  is  corroborated  by  brittle  failure  experiments  at  very  low  velocities 
[28].  Subsequently,  the  crack  propagation  takes  place  within  the  envelope  of  60°  —  70°,  which  is  in  agreement 
with  studies  in  [50,  51,  76].  The  dynamic  conditions,  as  well  as  boundary  constraints  are  responsible  for  the 
small  difference  between  X-VCFEM  results  and  those  in  [76].  Volume-averaged  or  macroscopic  shear  stress- 
shear  strain  behavior  for  this  problem  is  plotted  in  figure  4.15(b).  The  volume  averaging  of  the  local  stress 
and  strain  fields  over  the  entire  microscopic  domain  fl  is  performed  as 


CTy  (f ) 


tij  (t) 


(4.73) 


where  xj.  and  t  are  the  spatial  coordinates  and  time  (cumulative  increments  in  these  problems)  respectively, 
and 
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aij{t)  represents  the  effective  strain  field  caused  by  the  possible  displacement  jump  at  the  crack.  It  is 
calculated  along  the  crack  path  with  [ui{t)]  denoting  the  displacement  jump. 

Crack  propagation  in  sheared  plate  with  a  central  crack 

This  example  is  based  on  a  classical  problem  of  a  single  crack  propagation  in  a  large  plate  with  a  central  crack. 
The  plate  is  subjected  to  a  far  field  shear  load.  The  problem  was  experimentally  studied  by  Erdogan  and  Sih 
[28]  and  an  optical  micrograph  of  their  cracked  specimen  is  shown  in  figure  4.16(a).  The  specimen  material  in 
their  experiment  were  assumed  to  be  homogeneous,  isotropic  and  linearly  elastic  and  the  crack  was  assumed 
to  be  brittle.  A  single  element  of  dimension  10  m  x  8  m  in  X-VCFEM  is  used  to  simulate  this  experiment  as 
shown  in  figure  4.16(b).  The  initial  crack  length  is  a=1.6  m.  The  material  parameters  are:  Young’s  modulus 
E  =  100  GPa,  Poisson  ratio  v  =  0.3  and  the  cohesive  law  parameters  are:  (Tmax  =  0.1  MPa,  /?  =  1,  and 
de  =  1  X  lO^^m.  The  shear  stress  applied  on  the  top  and  bottom  surfaces,  is  varied  from  0  to  0.041  GPa  with 
plane  stress  assumptions.  As  shown  in  figure  4.16(b),  the  crack  path  predicted  by  X-VCEEM  compares  well 
with  the  observations  in  [28].  Figure  4.16(c)  shows  the  growth  of  the  crack  opening  displacement  components 
at  the  right  tip  A.  The  entire  computational  process  took  20  minutes  on  a  single  CPU  in  the  Pentium  4 
cluster  with  2.4Ghz  Intel  P4  Xeon  processors. 

Crack  propagation  in  three-point  bending  specimen 

Two  numerical  examples  are  considered  for  this  specimen.  In  the  first  example,  symmetric  mode  I  crack 
propagation  in  a  three-point  bending  test,  as  shown  in  figure  4.17,  is  modeled.  Plane  stress  conditions  are 
assumed  in  the  simulation.  This  problem  of  cohesive  crack  propagation  has  been  studied  by  Carpinteri  [14] 
using  node  release  technique  and  by  Mods  and  Belytschko  [59]  using  the  extended  FEM  or  XFEM.  The 
geometrical  dimensions  for  the  specimen  in  figure  4.17  are  b=  0.15  m,l=4b,  t(specimen  thickness)=  b,a=0, 
and  d=0.001  m.  The  material  properties  are:  Young’s  modulus  E  =  36,500  MPa,  Poisson  ratio  v  =  0.1, 
and  the  cohesive  parameters  are  a^ax  =  3.19  MPa,  and  /3  =  0.  The  X-VCFEM  solution  is  compared 
with  that  in  [59]  through  the  load-deflection  curve  of  figure  4.18.  The  cohesive  displacement  parameters  are 
6e  =  3.134796  x  10^®  m  and  6e  =  6.26959  x  10^®m  for  figures  4.18(a)  and  4.18(b)  respectively.  A  sharper 
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snap-back  is  seen  for  the  latter  case.  Excellent  match  is  observed  between  the  X-VCFEM  and  XEEM  results 
The  second  example  shows  a  mixed-mode  cohesive  crack  propagation  in  a  three-point  bend  test  due  to  an 
unsymmetrically  positioned  initial  crack.  The  problem,  shown  in  figure  4.19(a),  has  been  studied  by  Marian! 
and  Perego[56]  using  XEEM  under  plane  stress  conditions.  The  initial  crack  position  is  determined  by  the 
offset  ratio  a,  defined  as  the  ratio  of  the  distance  of  the  initial  crack  from  the  mid-span  cross-section  to  half 
of  the  beam  span.  The  material  Young’s  modulus  E  =  31370  MPa,  and  Poisson  ratio  v  =  0.2.  The  cohesive 
model  parameters  are  a^ax  =  4.4  MPa,  6e  =  0.07719298  mm  and  /3  =  1.0.  Once  again,  the  entire  domain 
is  represented  by  a  single  element  in  X-VCEFM  with  154  nodal  degrees  of  freedom.  Figure  4.19(b)  shows 
the  load-deflection  curve  for  two  values  of  the  offset  parameter,  i.e.  a  =  0.5  and  a  =  0.25.  The  initial  elastic 
response  in  the  load  P-displacement  u  curve  is  stiffer  and  also  the  peak  load  is  higher  for  higher  values  of 
a.  The  load-displacement  response  exhibits  softening  in  the  later  stages  of  crack  propagation  due  to  the 
significantly  evolved  crack.  The  path  of  crack  propagation  for  the  two  cases  are  shown  in  figures  4.19(c)  and 
(d).  The  cracks  move  towards  the  point  of  applied  load  and  align  themselves  perpendicular  to  the  edge  of 
the  specimen. Excellent  agreement  is  obtained  between  the  results  by  X-VCFEM  and  in  [56]. 

4.6.3  Mesh  independence  of  crack  propagation  with  X-VCFEM 

A  panel  with  domain  5  cm  x  3  cm  containing  two  initial  cracks  is  remotely  loaded  in  tension  as  shown  in 
figure  4.20(a).  The  problem  has  been  solved  by  Sharma  et.  al.  [87]  using  the  element  free  Galerkin  meshless 
method.  For  X-VCFEM  solution,  the  domain  is  meshed  into  two  elements  with  three  different  topologies 
shown  in  figure  4.20.  Plane  stress  conditions  are  again  assumed.  A  total  of  11  increments  is  used  to  model 
the  entire  crack  propagation  process.  The  material  parameters  are:  Young’s  modulus  E  =  207  GPa  and 
Poisson  ratio  n  =  0.3  and  cohesive  zone  parameters  are:  amax  =  0.1  MPa,  dg  =  1  x  10^®cm,  and  /3  =  1. 
The  three  figures  4.20(b,c,d)  show  no  mesh  dependence  of  the  X-VCFEM  predictions  and  the  comparison 
with  results  in  [87]  is  excellent. 
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4.6.4  Effect  of  cohesive  parameters  on  crack  evolution 


This  example  is  intended  to  investigate  the  effect  of  cohesive  parameters  on  crack  growth.  Cohesive  zone 
model  parameters,  e.g.  <Jmax  and  5^  in  equation  (4.11),  can  significantly  affect  the  propagation  and  overall 
behavior  of  a  cracking  material.  The  effects  of  these  cohesive  parameters  are  studied  for  crack  propagation 
in  a  sheared  plate  with  a  central  crack  subjected  to  a  far  field  shear  load.  This  classical  problem  was 
experimentally  studied  by  Erdogan  and  Sih  [28]  and  an  optical  micrograph  of  their  cracked  specimen  is  shown 
in  figure  4.21(a).  The  specimen  material  in  their  experiment  was  assumed  to  be  homogeneous,  isotropic  and 
linearly  elastic  and  the  crack  was  assumed  to  be  brittle.  A  single  element  of  dimension  10  m  x  8  m  in  X- 
VCFEM  is  used  to  simulate  this  experiment  as  shown  in  figure  4.21(b).  The  initial  crack  length  is  lo=l-6  m. 
The  material  parameters  are:  Young’s  modulus  E  =  100  GPa,  Poisson  ratio  v  =  0.3.  Eive  different  sets  of 
cohesive  parameters,  illustrated  in  figure  4.21(b)  are  considered  for  this  example.  These  are 

•  A:  CTtoox=3.0  MPa,  de=3.0  e-4m,  /3  =  1.0 

•  B:  CTtoox=6.0  MPa,  de=l-5  e-4m,  /3  =  1.0 

•  C:  (Ttoox=3.0  MPa,  de=6.0  e-4m,  (3  =  1.0 

•  D:  (Ttoox=6.0  MPa,  de=3.0  e-4m,  (3  =  1.0 

•  E:  (Ttoox=1-5  MPa,  de=6.0  e-4m,  (3  =  1.0 

As  shown  in  figure  4.21(b),  all  the  cases  correspond  to  the  same  cohesive  energy.  The  load  is  applied  by 
controlling  the  opening  of  crack  propagation  through  fixed  values  of  the  increment  A/  in  equation  (4.54). 
Further  a  uniform  shear  load  per  unit  length  tq  is  applied  on  the  top  and  bottom  surfaces  as  shown  in  figure 
4.21(c).  In  each  increment,  the  applied  load  is  scaled  by  the  arc-length  parameter  A  of  equation  (4.53),  to 
yield  an  equilibriated  applied  load  corresponding  to  a  prescribed  crack  propagation  length.  The  crack  path 
for  all  the  different  cohesive  parameters  predicted  by  X-VCFEM  are  very  similar  and  compare  well  with 
experimental  observations  in  [28].  However  a  considerable  dependence  on  cohesive  parameters  is  seen  in  the 
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shear-crack  length  response,  demonstrated  in  figure  4.21(d),  where  the  normalized  crack  length  is  defined  as 


The  current  crack  length 
The  initial  crack  length 


(4.75) 


This  points  to  the  fact  that  the  rate  of  propagation,  and  not  the  direction,  is  dependent  on  the  parameters 
for  this  problem.  For  the  cases  with  larger  peak  traction  cases:  B  and  D,  higher  applied  loads  are  needed 
for  causing  similar  crack  growths  as  for  cases  with  lower  peak  traction:  A  and  C.  Comparison  of  the  results 
for  cases  B  and  D,  show  that  a  smaller  dg  (case  B)  results  in  quicker  reduction  of  the  local  cohesive  traction. 
This  makes  the  overall  load  for  the  case  B  to  increase  slower  than  that  for  case  D  with  a  higher  dg.  The  case 
E  consistent  with  the  trends  exhibited  by  the  other  load  cases.  Although,  the  simulation  results  show  that 
both  Umax  and  dg  affect  the  crack  growth,  comparison  of  cases  A,  B,  C  with  D  shows  that  the  crack  growth 
is  more  sensitive  to  Umax  than  to  dg .  The  results  also  imply  that  the  cohesive  energy,  or  effectively  the  energy 
release  rate  Gc  does  not  alone  determine  the  properties  of  crack  propagation.  The  individual  parameters, 
affecting  the  shape  of  the  cohesive  law,  play  an  important  role  in  predicting  the  growth  characteristics.  These 
effects  are  also  tested  for  multiple  crack  growth  in  the  next  set  of  examples. 


4.6.5  Propagation  of  multiple  pre-existing  cracks 

The  final  set  of  examples  looks  into  the  growth  of  multiple  pre-existing  cracks  to  comprehend  the  effect  of 
morphology,  e.g.  distribution,  orientation  etc.. 

Firstly,  a  plate  with  five  randomly  located  cracks  is  simulated  under  a  tensile  loading  as  shown  in  figure 
4.22(a).  The  plate  has  dimensions  0.6  m  x  0.4  m;  material  parameters:  Young’s  modulus  E  =  10®  MPa 
and  Poisson  ratio  v  =  0.3;  and  cohesive  parameters:  Umax  =  0.1  MPa,  fi  =  1,  and  dg  =  1  x  10^®cm.  Figure 
4.22(b)  shows  the  final  positions  of  the  cracks  that  have  grown  with  the  loading.  The  cracks  propagate  across 
element  boundaries  and  are  attracted  to  each  other  in  certain  regions  till  they  nearly  merge. 

A  plate  with  28  randomly  located  and  oriented  cracks  is  simulated  under  a  tensile  loading.  Figures  4.23(a) 
and  (b)  show  the  two  microstructures  with  different  crack  distributions.  For  the  microstructure  1,  all  the 
cracks  of  equal  length  are  oriented  horizontally  and  their  distribution  is  random.  The  microstructure  2  has 
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cracks  of  random  length  and  orientation.  In  addition,  it  contains  a  cluster  of  8  cracks  in  a  otherwise  random 
distribution  as  shown  in  figure  4.23(b).  The  plate  is  of  dimension  0.1  m  x  0.1  m,  and  the  material  parameters 
are:  Young’s  modulus  E  =  10^  MPa  and  Poisson  ratio  p  =  0.3.  To  understand  the  effect  of  cohesive 
parameters  on  crack  propagation,  two  different  sets  of  cohesive  parameters  are  considered.  They  are: 

CP-1:  a^ax=E0  MPa,  de=l-0  e-5m,  /3  =  0.707 
CP-2:  a^ax=‘2.0  MPa,  de=0.5  e-5m,  /3  =  0.707 

A  uniform  tension  load  per  unit  length  a  is  applied  on  the  top  and  bottom  surfaces  as  shown  in  figure 
4.23(a,b).  In  each  increment,  the  applied  load  is  scaled  by  the  arc-length  parameter  A  of  equation  (4.53),  to 
yield  an  equilibriated  applied  load  corresponding  to  a  prescribed  crack  opening  deformation. 

Figures  4.24(a,b)  and  (c,d)  show  the  contour  plots  of  the  microstructural  stress  ciyy  together  with  evolved 
position  of  the  cracks  at  the  final  stage  of  loading,  for  the  two  sets  of  microstructures  and  cohesive  parameter 
respectively.  The  growth  pattern  of  each  crack  can  be  observed  by  comparing  with  its  initial  configuration 
in  figures  4.23(a)  and  (b).  The  cracks  propagate  across  element  boundaries,  interact  with  each  other  and  in 
some  cases,  they  merge.  The  relation  of  the  propagation  of  multiple  cracks  to  the  morphology  and  cohesive 
parameters  is  in  general  complicated.  However,  several  observations  can  be  made  based  on  the  results  of  the 
simulation  by  this  model. 

•  Larger  stress  concentrations  develop  at  tips  of  cracks  that  are  nearly  perpendicular  to  the  direction  of 
loading.  Consequently,  this  subset  of  cracks  grow  more  easily  than  others  that  are  more  aligned  with 
the  loading  direction.  From  figure  4.24(b)  and  (d),  it  can  be  seen  that  some  cracks  that  are  nearly 
parallel  to  the  load  direction  never  propagate. 

•  Stress  concentrations  are  higher  at  tips  of  longer  cracks.  The  reason  is  stress  concentrations  at  crack 
tips  come  from  the  external  load,  which  cannot  be  handled  by  the  weak  crack.  Longer  cracks  lead  to 
more  external  load  concentrating  to  tips.  This  is  verified  by  results  shown  in  figure  4.24(b)  and  (d), 
where  longer  cracks  are  easier  to  propagate  than  shorter  ones. 

•  Irrespective  of  the  initial  orientation,  the  evolved  crack  path  tends  to  align  in  a  direction  perpendicular 
to  the  applied  load  direction.  This  correspond  to  an  optimal  direction  for  releasing  the  cohesive  energy. 
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This  observation  is  dominant,  when  the  influence  of  nearby  cracks  on  the  local  stress  held  is  small. 
The  local  stress  held  for  this  phenomenon  is  mainly  governed  by  the  influence  of  applied  load  on  this 
single  crack. 

•  Cracks  are  attracted  towards  weak  surfaces,  such  as  other  cracks  or  voids  and  prefer  to  propagate  in 
those  directions.  This  may  be  attributed  to  the  fact  that  the  cohesive  energy  in  the  direction  of  these 
weaker  surfaces  with  lower  (or  zero)  tractions  is  naturally  lower  in  comparison  with  other  directions. 

•  Figures  4.24(b)  and  (d)  show  that  the  longest  crack  does  not  necessarily  evolve  from  a  cluster.  Not  all 
cracks  in  a  cluster  grow  considerably.  This  is  somewhat  in  contrast  to  observations  made  with  particle 
reinforced  composites,  where  almost  always  clusters  cause  a  local  stress  concentration.  The  interaction 
between  neighboring  cracks  contributes  to  the  enhancement  or  mitigation  of  stresses,  depending  on 
their  orientations  and  length.  This  dictates  their  propagation,  and  just  being  in  a  cluster  does  not 
guarantee  significant  growth. 

•  The  different  cohesive  parameters  show  very  little  difference  in  the  final  configuration  and  hence  the 
propagation  direction.  However,  the  rate  of  crack  growth  varies  considerably  with  these  parameters  as 
seen  in  the  crack  length-macroscopic  strain  plot  of  figure  4.25 

Figure  4.26  shows  the  macroscopic  stress-strain  response  for  the  two  microstructures  and  cohesive  pa¬ 
rameters.  Even  before  the  cracks  propagate  (corresponding  to  the  change  in  slope),  the  stiffness  of  the 
microstructure  2  is  higher  than  that  of  microstructure  1  due  to  a  higher  level  of  effective  damage  caused  by 
crack  lengths  and  more  importantly  orientations.  Orientations  perpendicular  to  the  load  direction  causes  a 
larger  reduction  in  stiffness  in  comparison  with  other  directions.  With  additional  loading,  the  overall  damage 
caused  by  the  growth  of  cracks  is  also  higher  for  the  microstructure  1.  This  is  seen  by  the  lower  values  of  the 
macroscopic  stress  for  this  case.  The  effect  of  the  cohesive  parameters  on  the  stress-strain  response  is  quite 
pronounced.  The  maximum  macroscopic  stress  for  both  microstructures  increases  significantly  for  higher 
values  of  Umax^  even  though  the  cohesive  energy  is  the  same  for  the  two  cohesive  models.  This  is  caused  by 
a  slowdown  in  the  growth  rate  of  the  cracks  with  overall  deformation. 
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4.7  Concluding  Remarks 


The  extended  Voronoi  cell  finite  element  model  is  developed  in  this  chapter  to  predict  initiation  and  growth 
of  damage  by  crack  propagation  in  brittle  matrix.  The  cracks  are  modeled  by  a  linear  cohesive  zone  model 
and  their  incremental  directions  and  growth  lengths  are  determined  in  terms  of  the  cohesive  energy  near  the 
crack  tip.  Important  enhancements  are  made  to  the  element  to  allow  stress  discontinuities  across  the  cohesive 
crack  and  to  accurately  depict  the  crack  tip  stress  concentrations.  These  features  are  accommodated  through 
the  incorporation  of  (a)  branch  functions  in  conjunction  with  level  set  methods  across  crack  contours,  and 
(b)  adaptive  multi-resolution  wavelet  functions  in  the  vicinity  of  the  crack  tip.  Several  problems  are  solved 
and  compared  with  existing  solutions  in  the  literature  for  validation  of  the  X-VCFEM  algorithms,  both  with 
respect  to  macroscopic  (load-deformation  behavior)  and  microscopic  (crack  path).  The  X-VCFEM  results 
show  excellent  accuracy  in  their  comparison  with  analytical  and  other  numerical  solutions.  Also  compari¬ 
son  with  ABAQUS  shows  the  efficiency  of  X-VCFEM.  Numerical  simulations  are  conducted  with  different 
and  de  to  understand  the  effect  of  cohesive  parameters  on  the  crack  propagation.  It’s  observed  that  in 
addition  to  the  total  cohesive  energy,  the  individual  parameters  have  effects  on  crack  growth.  The  effect  of 
geometrical  information  of  multiple  pre-existing  cracks,  including  the  lengths,  positions  and  orientations  of 
cracks,  on  their  propagation  is  studied  by  simulating  a  plate  with  28  randomly  located  and  oriented  cracks. 
Simulation  results  show  that  the  crack  with  a  longer  length  and  nearly  perpendicular  to  load  direction  is 
easier  to  propagation  than  other  cracks.  Cracks  propagation  direction  is  dependent  on  the  local  stress  field, 
which  is  managed  by  both  the  external  load  and  nearby  material  phases,  such  as  other  cracks  in  a  cluster. 
This  research  reveals  the  significance  of  analyzing  large  regions  of  the  microstructure  and  proves  the  effec¬ 
tiveness  of  the  X-VCFEM.  The  simulation  results  based  on  X-VCFEM  could  also  provide  positive  feedback 
for  design  modification. 

Based  on  the  study  on  interfacial  debonding  and  cohesive  matrix  cracking  in  composites,  the  interaction 
of  the  two  damage  phenomena  is  studied  in  the  next  chapter,  where  the  X-VCFEM  is  improved  and  a 
criterion  for  assessing  the  direction  of  damage  development  is  proposed. 
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tra 

0.6 

0.6 

0.6 

0.5 

0.5 

0.5 

0.4 

0.4 

0.4 

h 

0.05 

0.1 

0.2 

0.05 

0.1 

0.2 

0.05 

0.1 

0.2 

wmmm 

1.019 

1.014 

1.040 

1.017 

1.014 

1.027 

1.020 

1.020 

1.035 

Table  4.1:  Normalized  stress  intensity  factors  {Kj/Kref)  for  different  values  of  tra  and  bi  in  the  multi¬ 
resolution  wavelet  representation. 


Case  1 

Case  2 

Case  3 

Case  4 

Case  5 

WSHBM 

0 

1 

2 

3 

4 

KilKref 

1.1642 

1.0361 

1.0208 

1.0062 

1.0020 

ASE 

96.45% 

45.91% 

7.06% 

3.01% 

Table  4.2:  Errors  with  varying  enrichment  order  of  multi-resolution  wavelet  functions  for  the  different  cases. 
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(a) 


Itm 


Figure  4.1:  (a)  A  mesh  of  Voronoi  cell  elements,  each  containing  a  single  pre-existing  crack,  (b)  a  typical 
Voronoi  cell  element  showing  different  topological  features  and  loads. 


(a)  (b) 

Figure  4.2:  Normal  and  tangential  traction-separation  behavior  for  the  linear  cohesive  zone  model. 
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Crack  tip  2 


Crack  Surface 


Crack  tip  1 


Figure  4.3:  (a)  A  schematic  diagram  of  a  crack  surface  showing  parameters  related  to  the  distance  functions; 
(b)  depiction  of  the  branched  stress  function  near  a  crack  for  s  =  0,  t  =  0. 
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VCFEM 


Regions  covered  by  wavelet 
bj*aj  basis  functions 
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Figure  4.4:  Distribution  of  multi-resolution  wavelet  bases  around  a  crack  tip:  (a)  Crosses  (x)  refer  to  the 
location  of  the  origin  of  the  basis  vectors  at  a  lower  level  corresponding  to  dilation  parameters  (fro  and  tVc) 
and  (b)  adaptively  upgraded  to  higher  level  wavelet  bases  with  the  addition  of  the  next  level  of  bases  at 
locations  indicated  by  the  (□). 
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No 


End 


Figure  4.6:  A  flowchart  of  the  solution  method. 
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Figure  4.7:  The  stress  at  x  =  —0.3  for  a  double  cantilever  beam  to  demonstrate  the  effect  of  the 
branched  stress  function. 


Figure  4.8:  Algorithms  for  incremental  propagation  of  cohesive  cracks:  (a)  for  direction  and  incremental 
length,  (b)  a  cohesive  crack  going  through  the  inter-element  boundary,  (c)  for  merger  with  other  cracks  and 
(d)  for  evaluation  of  J—  integral. 


84 


Figure  4.9:  Subdivision  of  the  Voronoi  cell  element  for  Gaussian  quadrature,  with  a  higher  density  of 
integration  points  near  the  crack  tip. 
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(a)  (b) 


Figure  4.10:  (a)  A  center  cracked  plate  loaded  in  tension,  (b)  a  single  X-VCFEM  element  with  prescribed 
boundary  conditions 
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Figure  4.12:  X-VCFEM  generated  stress  Oyy  at  y  =  0  for  different  enrichment  orders  of  the  wavelet  basis 
functions. 


88 


Remote  load  g  (MPa) 


(b) 


Figure  4.14:  (a)  A  plate  with  an  edge  crack  under  remote  tension  load,  (b)  comparison  of  load-deformation 
curves  by  X-VCFEM  and  ABAQUS. 


Macroscopic  shear  stress  a  (GPa) 


ra.i 
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Figure  4.15:  (a)  Prediction  of  the  crack  path  by  X-VCFEM  for  the  Kalthoff  experiment,  (b)  the  macroscopic 
stress-strain  response. 
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Figure  4.17:  A  three-point  symmetric  bending  specimen. 
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Figure  4.18:  Comparison  of  normalized  load-deflection  curves  for  the  three-point  bending  beam:  (a)  6, 
3.134796  X  10-®  m  and  (b)  de  =  6.26959  x  lO^®  m. 
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Figure  4.19:  (a)  A  three-point  bending  specimen  with  an  unsymmetric  initial  crack,  (b)  comparison  of 
load-deflection  curves  from  X-VCFEM  and  literature  [56]  ,  (c)  and  (d)  comparison  of  the  crack  paths  by 
X-VCFEM  with  that  in  [56]  for  a  =  0.25  and  a  =  0.5,  respectively. 
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(a)  (b) 


(c) 


(d) 


Figure  4.20:  A  plate  with  two  cracks  in  arbitrary  locations  modeled  by  X-VCFEM  using  elements  of  different 
topologies  located  cracks,  (b,c  and  d)  show  crack  path  at  the  end  of  the  loading  for  the  different  elements 
and  also  a  comparison  with  [87] . 
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(c) 


(d) 


Figure  4.21:  (a)  Optical  micrograph  showing  the  path  of  cracking  in  a  plate  with  a  central  crack  subjected  to 
far-held  shear  [28],  (b)  5  different  sets  of  cohesive  parameters  for  X-VCFEM  simulations,  (c)  corresponding 
crack  path  generated  by  X-VCFEM,  (d)  comparison  of  the  growth  of  cracks  for  different  cases. 
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Figure  4.23:  Crack  propagation  in  two  square  regions  containing  28  cracks  by  X-VCFEM:  (a)  domain  with 
horizontal  cracks  of  equal  length  and  random  distribution,  (b)  domain  with  random  orientation,  length  and 
distribution  of  cracks  but  containing  a  cluster. 
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Figure  4.24:  Crack  propagation  in  two  square  regions  containing  28  cracks  by  X-VCFEM:  (a,b)  contour  plots 
of  (jyy  (MPa)  with  cohesive  parameters  CP-1  for  the  domains  in  figure  23  (a)  and  (b),  (c,d)  contour  plots  of 
Oyy  (MPa)  with  cohesive  parameters  CP-2  for  the  domains  in  figure  23  (a)  and  (b) 
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Figure  4.26:  Macroscopic  stress-strain  response  for  different  microstructural  morphologies  and  cohesive  pa¬ 
rameters. 
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Chapter  5 


Extended  Voronoi  Cell  Finite  Element 
for  Modeling  Interfacial  Debonding 
with  Matrix  Cohesive  Cracking  in 
Fiber  Reinforced  Composites 

5.1  Introduction 

Interfacial  debonding  and  cohesive  cracks  propagation  in  brittle  matrix  are  two  important  damage  phenomena 
in  fiber-matrix  composites.  Experiments  show  that  the  two  damage  phenomena  appear  in  the  same  material, 
where  the  failure  often  starts  from  the  interface  between  fiber  and  matrix,  and  is  subsequently  advanced  into 
matrix.  Researches  regarding  a  crack  meeting  a  bimaterial  interface  to  either  deflect  along  the  interface  or 
penetrate  into  the  next  layer  were  made  in  [2,  41,  42,  55],  where  the  criterion  of  deflection  versus  penetration 
was  established  based  on  the  energy  release  rate  and  fracture  energy.  However,  the  present  research  is  aimed 
at  only  elastic  cases,  which  requires  that  the  fracture  process  zone  at  the  crack  tip  is  small  compared  to 
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the  size  of  the  crack  and  the  size  of  the  specimen.  In  chapter  4,  cohesive  zone  models  are  introduced  into 
VCFEM  to  study  damage  of  interface  and  matrix,  where  the  stress  field  in  composites  are  described  by  a 
set  of  specifical  functions  accurately.  And  the  effect  of  cohesive  parameters  and  morphological  distributions 
are  studied  as  important  factors  to  the  damage  process.  All  work  in  the  two  chapters  are  theoretical  and 
computational  preparation  for  solving  interfacial  debonding  problems  coupled  with  matrix  cracking. 

In  this  chapter,  a  criterion  based  on  cohesive  zone  models  is  proposed  for  assessing  the  direction  of  damage 
development.  The  improved  X-VCFEM  is  developed  for  modeling  both  the  growth  of  interfacial  debonding 
and  the  propagation  of  multiple  cohesive  cracks  in  the  brittle  matrix  of  fiber-reinforced  composites.  The 
mechanics  theories  and  numerical  algorithm  in  previous  chapters  are  organized  as  an  organic  whole,  not  just 
a  simple  superposition.  It  begins  with  the  X-VCEEM  formulation  and  numerical  implementation,  followed 
by  the  numerical  example  showing  the  effectiveness  of  this  model  and  the  interaction  of  interface  and  crack 
propagation. 

5.2  Extended  Voronoi  cell  FEM  formulation  for  composites  with 
interfacial  debonding  and  matrix  cracking 

The  Voronoi  cell  finite  element  mesh  for  a  microstructure  with  both  debonded  interfaces  and  cohesive  cracks 
is  shown  in  figure  5.1(a),  where  the  region  is  divided  into  an  unstructured  finite  element  mesh  of  arbitrary 
Voronoi  cells.  A  typical  Voronoi  cell  element  fie  is  shown  in  figure  5.1  (b).  Each  VC  element  is  composed 
of  the  matrix  phase  (fi™),  the  inclusion  phase  (fie),  the  interface  (fim),  and  cracks  (ficr),  such  that  fie  = 
fim  U  U  ^in  U  ^cr,  where  interface  and  cracks  are  consider  as  zero  thickness  regions.  The  element  outer 
boundary  consists  of  the  prescribed  displacement  boundary  (Tum),  prescribed  traction  boundary  (ri™)  and 
the  inter-element  boundary  (Eto),  so  i.e.  9fie  =  Tum  U  U  Compatible  displacement  conditions  apply 
on  9fie.  9fi°  has  an  outward  normal  n°  (=n™),  while  n®  is  the  outward  normal  to  Sfig.  In  order  to  describe 
debonding  with  progressing  deformation  through  decohesion,  the  interface  is  lined  with  a  set  of  node-pairs 
with  nodes  belonging  to  the  matrix  interface  (9fi™)  and  inclusion  interface  {d^%)  respectively.  The  traction 
fcoft  between  node-pairs  on  the  crack  surface  are  modeled  by  the  cohesive  zone  traction-separation  law.  The 
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behaviors  of  cohesive  cracks  in  the  brittle  matrix  are  described  by  a  similar  method,  where  nodes  in  node-pairs 

1  2 

are  arranged  at  different  sides  of  a  crack  (Fcr  and  Fcr).  In  the  incremental  assumed  stress  hybrid  X-VCFEM 
formulation,  the  complementary  energy  functional  for  each  element  is  expressed  in  terms  of  increments  of 
stress  and  displacement  fields  as: 


ne((T,  A(T,U,  Au)  = 
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-b 

-b 
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(5.1) 


Here  B  is  the  complementary  energy  density  and  the  superscripts  m  and  c  correspond  to  variables  asso¬ 
ciated  with  the  matrix  and  inclusion  phases,  cr™  and  <7°  are  the  equilibrated  stress  fields,  e™  and  e°  the 

1  2 

corresponding  strain  fields  in  different  phases  of  each  Voronoi  element.  Also,  u®,  u™,  u®,  u®*"  and  u®*"  are 

1  2 

the  kinematically  admissible  displacement  fields  on  d^lg,  9fl™,  Ter  and  Ter  respectively.  The  prefix  A 

corresponds  to  increments.  The  term  in  the  box  in  equation  (5.1)  provide  the  work  done  by  the  interfacial 

tractions  T™  =  T™n™  -b  T™t™  due  to  interfacial  separation  (u™  —  u®),  where  T™  and  T™  are  the  normal 

and  tangential  components  that  are  described  by  cohesive  laws  at  the  interface.  Similarly,  the  last  term 

provide  the  work  done  by  the  cohesive  tractions  T®*"  =  T®’'n®’'  -b  Tj®’’!®’’  due  to  displacement  separation 
1  2 

(u®*"  —  u®*")  along  the  crack,  where  and  Tj®’’  are  the  normal  and  tangential  components  of  the  cohesive 
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force.  The  total  energy  for  the  entire  composite  domain  is  obtained  by  adding  the  energy  functionals  for  N 


elements  as 


JV 

n  =  Ene 

e=l 


(5.2) 


5.2.1  General  element  assumptions  and  weak  form 

In  the  absence  of  body  forces,  two  dimensional  stress  fields  satisfying  equilibrium  relations  can  be  generated 
from  the  Airy’s  stress  function  $(x,y).  In  the  incremental  formulation,  stress  increments  in  matrix  and 
inclusion  are  obtained  from  derivatives  of  the  stress  functions  A$™(x,y)  and  A^^{x,y)  as: 
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=  [P™(x,y)]{A/3™}, 


=  [P%x,y)]{An 


(5.3) 


where  {A/3™}  and  {A^^}  are  the  column  of  unknown  stress  increment  coefficients.  Convergence  properties 


and  efficiency  of  X-VCFEM  depend  on  the  choice  of  $™.  These  functions  should  adequately  account  for  the 


geometry  and  location  of  the  heterogeneity  in  the  element,  so  stress  functions  for  matrix  are  decomposed 


into  (a)  a  purely  polynomial  function  (b)  a  reciprocal  function  $™j.!  (c)  a  branch  function  ^^anch 

and  (d)  wavelet  functions  $™^„  (^™  =  $™;^  +  $™,  +  ^Zanch  +  The  selection  of  stress  functions 


are  discussed  in  chapter  4  detailed.  Inclusion  stress  functions  are  admitted  as  polynomial  function 
(^c  _  Compatible  displacement  fields  satisfying  inter-element  continuity  on  the  element  boundary 

and  intra-element  continuity  on  both  the  interface  and  the  crack  face  Ter  are  generated  by 
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interpolation  of  nodal  displacements  as: 


{Au®}  =  [L®]{Aq®}  onSOe 
{Au™}  =  [L™]{Aq™}  on  50™ 

{Au"}  =  [L=]{Aq=}  on  50^ 

{A^’-}  =  [L'’-]{Aq“’’}  on  r',, 

=  [L'niAq"’’}  on  (5.4) 

1  2 

The  interpolation  matrices  [L®],  [L™],  [L®],  [L®’’],  [L®’’]  for  the  nodal  displacements  on  the  respective  bound¬ 
aries  are  constructed  using  standard  linear  or  hierarchical  shape  functions.  Since  nodes  on  the  inter¬ 
face  and  crack  surfaces  are  always  belonging  to  some  node-pair,  the  interpolation  matrices  are  chosen  as 
[L™]  =  [L®]  and  [L®’’]  =  [L®’’]. 

Substituting  the  interpolations  of  stress  and  displacement  fields  from  equations  (5.3)  and  (5.4)  into  equation 
(5.1)  results  in  the  matrix  form  of  the  element  complimentary  energy 
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(5.5) 
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where 


[H™]  = 

[  [P™]^[S™][P™]dO™  [H®]  = 

[  [P®]"’[S®][P®]dO, 

'fie 

[G®]  = 

[  [P™]^[n®][L®]d50™  [G™]  = 

[  [P™]^[n™][L™]d5a 

Jon^ 

Jan^ 

[G®]  = 

[  [P®]^[n®][L®]d50c  [G"1  = 

L  [P]^[n®’’][L,,]dr,, 

Jan^  J 

r^r- 

[G®’’]  = 

L  [P]^[n®’’][Ll]dr,,  {t}  =  f 

■frer  JTtrr, 

[Le]^{t  -I-  AtjdrtTO 

Construction  of  appropriate  stress  functions  with  optimally  high  resolution  is  necessary  for  accurately  de¬ 
picting  high  stress  gradients  near  the  crack  tip. 


5.2.2  Solution  Method 


Crack  growth  in  multiply  cracked  materials  is  solved  using  an  incremental  approach,  where  a  set  of  elemental 
and  global  equations  are  solved  in  each  increment  for  stresses  and  displacements. 

1.  Local  equations  for  each  element  are  obtained  by  setting  the  variation  of  equation  (5.5)  with  respect 
to  the  stress  coefficients  A/3'^  and  A/3‘^  to  zero.  This  results  in  the  weak  form  of  the  element  kinematic 
relations 
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or  in  a  condensed  form 


[H]e{;a  +  A^}e  =  [G]e{q  +  Aq}e 


(5.7) 


(5.8) 
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Since  equation  (5.8)  is  linear,  the  stress  coefEcients  can  be  directly  expressed  in  terms  of  the  nodal  displace¬ 
ments,  provided  the  element  [H]e  matrix  is  invertible. 

2.  Subsequently,  the  weak  forms  of  the  global  traction  continuity  conditions  are  solved  by  setting  the  varia¬ 
tion  of  the  total  domain  energy  functional  with  respect  to  the  generalized  displacement  components  to  zero. 
This  results  in  the  weak  form  of  the  traction  reciprocity  conditions 


t 


1  2 

T 

N 

E 

[Qe]  [Qj  [Qcrj  _[Qcrj 

II 

e=l 

_  [0]  [0]  [G=]  [0]  [0] 

e 

<  J 

e=l 

e 

f*cr 

^coh 


(5.9) 


or  in  a  condensed  form: 


^[G]n/3  +  A/3}e=£{Te.i}e 
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(5.10) 


The  forces  at  the  interface  and  crack  surface  are  expressed  in  terms  based  on  the  cohesive  energy 


as 
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Combining  equations  (5.8)  and  (5.10)  and  eliminating  the  stress  coefficients  {y5+Ay3}e,  results  in  the  equation 
for  solving  the  generalized  displacements 


N 


N 


^{[G]nH]-nG]e}{q  + Aq}  = 


(5.13) 


Equation  (5.13)  is  a  nonlinear  matrix  equation  system  due  to  the  cohesive  laws.  Consequently,  a  Newton- 
Raphson  iterative  solver  is  invoked  to  solve  for  the  increments  of  nodal  displacements.  The  linearized  form 
of  equation  (5.13)  for  the  j-th  iteration  is 

j -  E{[G]^[H]-MG]}e{q  +  Aq}|  (5.14) 

U=1  e=l  J 


which,  in  a  condensed  form  is 


[K^Ydq^  =  {RLJ  -  {RL}^'  (5.15) 

Many  numerical  examples  in  Chapter  4  prove  that  only  a  Newton-Raphson  iterative  solver  cannot  obtain 
the  entire  failure  solution  for  the  problems  with  damage,  especially  when  a  snap-back  appears  in  the  load- 
deformation  curve. 

According  to  the  arc-length  method  proposed  in  [19,  20,  84],  an  unknown  loading  parameter  (A  -I-  dX)  is 
introduced  to  govern  the  load  increments.  Equation  (5.15)  is  modified  with  this  loading  parameter  as 

[K^YdcY  =  {X^  +  dA^HRLJ  -  {RltV  (5.16) 

where  both  dX^  and  dq-^  are  unknowns,  and  dX^  can  be  either  positive  or  negative.  The  orthogonality 
condition  is  chosen  to  be  the  constraint  equation  required  by  the  additional  unknown  dX^ . 


no 


5.2.3  Stability  conditions 


Following  the  stability  conditions  derived  for  displacement-based  and  stress-based  finite  element  approxima¬ 
tions  in  [4,  12,  104],  the  stability  conditions  of  the  stress-displacement  field  variational  problem  in  X-VCFEM 
are  stated  in  section  4.2.3.  They  are  positive  definite  [H™]  and  m,  unique  stress  interpolation  functions, 
and  non-zero  stress  parameters  for  all  non-rigid  body  displacement  fields.  The  two  conditions  can  be  sat¬ 
isfied  by  implementing  numerical  methods  in  section  4.5.  And  the  third  one  is  accomplished  by  choosing 
>  n®  -I-  n™  -I-  n®’’  *  2  -  3  and  npc  >  n®  -  3. 

5.3  Aspects  of  Numerical  Implementation 

5.3.1  Adaptive  criteria  for  cohesive  crack  growth 

A.  The  criterion  for  the  incremental  cohesive  crack  advance  into  matrix: 

The  static  deflection/penetration  behavior  at  an  interface  has  been  the  subject  of  numerical  research  efforts 
in  the  past  years  and  many  significant  results  for  various  kinds  of  materials  have  been  obtained([2,  41, 
42,  55]).  The  fracture  toughness  ratio  of  the  interface  and  the  matrix  material  has  been  identified  as  the 
most  important  parameter  governing  the  crack  deflection/penetration  phenomenon.  Predicting  crack  growth 
requires  to  calculate  the  energy  release  rate,  G,  and  a  knowledge  of  the  surface  fracture  energy,  Gc-  In  this 
chapter,  we  denote  by  G*  and  G*  the  energy  release  rate  and  the  critical  energy  release  rate  for  the  case  of 
growth  along  interfaces,  and  by  G™  and  G™  the  corresponding  quantities  for  penetration  into  matrix. 

As  seen  in  previous  results,  stress  concentration  always  appears  in  the  matrices  around  fibers,  which  results 
in  cohesive  interfaces  between  fiber  and  matrix  becoming  weak  and  even  debonded.  Simultaneously,  damage 
at  the  interface  results  in  larger  concentrated  stress  fields  in  matrix.  Once  the  stress  in  matrix  reaches  some 
critical  value,  the  material  at  this  matrix  point  might  become  softening  and  damage  propagates  into  matrix 
from  the  interface.  All  points  with  critical  stresses  are  regarded  as  the  candidate  damage  position,  where  the 
criterion  is  necessary  for  selecting  the  crack  growth  direction,  along  the  interface  or  branching  into  matrix. 
The  candidate  positions  are  usually  chosen  from  the  Gaussian  integration  points  on  the  interface.  In  the 
program,  18  Gaussian  integration  points  are  distributed  between  any  two  consecutive  nodes  at  the  interface. 


Ill 


At  the  candidate  points,  the  criterion  for  assessing  the  crack  penetrating  into  matrix  is  defined  as 


(5.17) 

In  this  thesis,  the  energy  release  rate  is  calculated  based  on  cohesive  zone  models,  which  are  shown  in  figure 
5.2.  The  bilinear  model  in  figure  5.2  (a)  is  for  describing  the  damage  at  interface,  and  the  linear  model  in 
figure  5.2  (b)  is  for  the  matrix  cracking.  The  areas  of  the  shadow  regions  express  the  current  energy  release 
rates  G*  and  G™.  According  to  the  relation  between  the  cohesive  energy  (j)  for  complete  decohesion  and  the 
critical  energy  release  rate  Gc  in  equation  (4.55),  the  critical  release  rates  for  interface  and  matrix  are 

Gi  =  \ol,Ji  and  G^  =  ^aZaJT  (5-18) 

In  order  to  obtain  the  energy  release  rate  G™,  the  effective  cohesive  traction  t  is  calculated  according  to 
stresses  (ct^x,  (^yy  and  a^y)  at  every  candidate  point.  Recalling  equations  (4.58-4.62)  in  Chapter  4,  effective 
cohesive  traction  t{ac),  the  cohesive  energy  ^(ctc)  and  the  energy  release  rate  G™  are  obtained 

t{a)  = 

\  a  —  CTa,3/Sin2a  +  ayycos^a)^  +  ^CTa,xSin2a  +  axyCOs2a  +  ^ayySm2ay 

V  A  A 

(5.19) 


G™  = 

Xm 

=  ^(ac)  =  2  ®  {oZax  t{acf) 

^^max 

(5.20) 

where  is  the  angle  maximizing  the  cohesive  energy. 

The  current  energy  release  rates  for  interface,  G*,  is  obtained 

G*  =  < 

d  <  dc 

(5.21) 

2  (5^ -(5c  f 

5  >  5c 
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According  to  equations  (5.18,  5.20,  5.21)  and  inequality  (5.17),  the  damage  propagation  directions  are 
determined  at  the  candidate  points. 

B.  Direction  and  length  of  the  incremental  cohesive  crack  advance: 

Recalling  results  in  chapter  4,  the  direction  of  matrix  cracking  etc  is  obtained  at  the  damage  onset  points  as 


OLc  — 


arctan 


20-* 


(5.22) 


The  sign  in  equation  (5.22)  is  chosen  as  the  one  that  maximizes  the  cohesive  energy  by  satisfying  the 
condition  in  equation  (4.61)  c. 

Upon  establishing  the  direction  of  incremental  cohesive  crack  growth  etc,  the  length  of  cohesive  zone  advance 
(A/)  should  be  estimated  in  the  crack  evolution  scheme  according  to  the  same  algorithm  shown  in  chapter 
4  as: 


Al=  \AB\,  (5.23) 

9  A  —  9b 

where  B  is  a  point  close  to  A  in  the  direction  of  crack  propagation. 

5.3.2  Generation  of  [G^] 

Once  damage  is  driven  from  interface  into  matrix,  two  node-pairs  (mi,  ni)  and  (m2,  712),  shown  in  figure 
5.3,  are  added  at  the  interface,  where  nodes  mi  and  m2  are  at  the  matrix  side  and  nodes  rii  and  712  are 
at  the  inclusion  side.  The  separation  between  mi  and  m2  describes  the  displacement  discontinuity  at  the 
crack  surface.  Since  crack  doesn’t  propagates  into  the  inclusion,  the  node-pair  (rii,  712)  merges  by  sharing 
the  same  displacement.  This  can  be  implemented  at  assembling  matrix  [G°].  In  matrix  [G^’],  the  elements  in 
column  DOFn2  are  added  to  the  corresponding  elements  in  column  DOFni,  and  the  entire  column  DOFn2 
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is  assigned  zero.  The  process  is  shown  in  equation  (5.24)  as 


DOF^ 


V- 


/ 


DOFr, 


*  +  # 


••  *  +  # 

V---  *  +  # 


D0F„2 

# 

# 

# 

D0F^2 

0 

0 

0 


=> 


n^c  XTlqC 


n^c  XTlqC 


5.4  Numerical  Example 


(5.24) 


An  example  with  a  square  microstructure  containing  a  single  circular  fiber  with  a  debonding  interface  is 
considered  to  check  the  effectiveness  of  X-VCFEM  and  study  the  interaction  between  the  interface  and 
matrix  cohesive  cracking.  The  geometrical  dimensions  for  the  specimen  in  figure  5.4(a)  are  a=  20  mm, 
r=  5  mm.  The  material  parameters  for  matrix  and  fiber  are:  Young’s  modulus  =  72  GPa,  Ef  = 
450  GPa,  and  Poisson  ratio  Vm  =  0.32,  Vf  =  0.17,  where  subscript  (•)„  and  (•)/  denote  matrix  and 
fiber  respectively.  The  interface  uses  the  bilinear  cohesive  zone  model  with  the  properties  =  0.04  GPa, 
6c  =  0.001mm,  d*  =  0.02mm,  /3*  =  0.707.  The  linear  cohesive  zone  model  is  used  to  describe  matrix  cracking 
with  parameters:  =  0.05  GPa,  d™  =  0.002mm,  and  /3™  =  1.  The  cohesive  parameters  are  chosen 

to  make  so  that  the  damage  starts  from  interface  instead  of  matrix.  Under  plane  strain 

conditions,  the  displacement  boundary  conditions  are  shown  in  figure  5.4(a).  The  whole  microstructure  is 
modeled  with  one  X-VCFEM  element,  consisting  of  16  nodes  on  the  cell  boundary  and  20  node  pairs  on  the 
interface  for  displacement  interpolation.  Before  damage  propagates  into  matrix,  the  stress  functions  in  this 
example  consist  of  102  terms  of  polynomial  functions  and  45  terms  of  reciprocal  functions.  After  the  cracks 
advance  into  matrix,  one  branch  function  and  16  wavelet  functions  are  added  into  the  stress  interpolation  for 
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each  crack.  Figures  5.4(b)  shows  the  contour  plots  of  the  microstructural  stress  ayy  together  with  evolved 
position  of  the  cracks  at  the  final  stage  of  loading.  The  growth  pattern  of  each  crack  can  be  observed  by 
comparing  with  its  initial  configuration  in  figures  5.4(a),  where  there  is  no  matrix  cracks. 

The  relation  of  the  propagation  of  multiple  cracks  to  the  interface  debonding  is  in  general  complicated. 
However,  several  observations  can  be  made  based  on  the  results  of  the  simulation  by  this  model. 

•  As  shown  in  figures  5.4(b),  the  lower  stress  at  point  A  implies  that  interface  there  becomes  weak  even 
debonded,  which  results  in  the  load  can  not  be  transfered  into  the  fiber  at  this  position  effectively.  The 
positions  with  concentrated  stress  bifurcate  from  point  A  and  move  to  left  and  right  sides  respectively 
along  the  interface,  which  might  drive  the  damage  into  matrix  from  the  interface.  The  same  thing 
happens  at  the  bottom  point  of  the  circular  fiber. 

•  Due  to  symmetry,  four  cracks  propagate  into  matrix  from  the  interface.  Largest  stresses  appear  at  tips 
of  cracks  and  the  stress  in  fiber  is  released.  In  this  example,  since  cracks  result  in  larger  concentrated 
stress  than  interfaces,  cracks  propagation  in  matrix  becomes  the  key  damage  phenomenon  in  following 
failure  process. 

•  The  evolved  crack  path  tends  to  align  in  a  direction  perpendicular  to  the  applied  load  direction,  which 
agrees  with  the  observation  in  chapter  4. 

5.5  Concluding  Remarks 

The  extended  Voronoi  cell  finite  element  model  is  improved  in  this  chapter  to  predict  the  damage  advanc¬ 
ing  into  matrix  and  study  the  interaction  between  interfacial  debonding  and  matrix  cracking.  Polynomial 
functions,  reciprocal  functions,  branch  functions  and  wavelet  functions  are  made  to  the  element  stress  in¬ 
terpolations  to  accurately  depict  the  stress  discontinuities  and  concentrations  at  interfaces  and  cracks.  The 
damage  in  interface  and  matrix  are  modeled  by  cohesive  models.  A  criterion  for  assessing  the  crack  pen¬ 
etrating  into  matrix  is  proposed,  which  is  based  on  the  energy  release  rate  and  cohesive  energy.  A  square 
specimen  containing  a  single  circular  fiber  with  a  debonding  interface  is  considered  to  check  the  effectiveness 
of  X-VCFEM  and  the  criterion. 
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The  damage  analysis  in  fiber-reinforced  composites  is  in  general  complicated.  X-VCFEM  is  easy  to  be  ex¬ 
tended  to  study  effects  of  material  properties  and  geometric  characterization,  such  as  clustering,  alignment, 
fiber  shape,  relative  sizes  etc.,  which  are  critical  to  the  failure  process  in  the  microstructure.  This  will  be 
explored  in  the  future  work. 


116 


Chapter  6 


Concurrent  Multi-level  Model  for 
Damage  Evolution  in 
Microstructurally  Debonding 
Composites 

6.1  Introduction 

Analysis  of  composite  materials  with  microstructural  heterogeneities  is  conventionally  done  with  macroscopic 
properties  obtained  by  homogenizing  response  functions  in  the  representative  volume  element  (RVE)  from 
microscopic  analyses  at  smaller  length  scales.  While  these  “bottom-up”  homogenization  models  are  efficient 
and  can  reasonably  predict  macroscopic  or  averaged  behavior,  such  as  stiffness  or  strength,  they  have  limited 
predictive  capabilities  with  problems  involving  localization,  failure  or  instability.  Assumptions  of  macroscopic 
uniformity  and  RVE  periodicity,  the  two  basic  requirements  of  homogenization,  break  down  under  these 
circumstances.  The  uniformity  assumption  ceases  to  hold  in  critical  regions  of  high  local  solution  gradients. 
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such  as  near  free  edges,  interfaces,  material  discontinuities  or  evolving  damage.  RVE  periodicity,  on  the  other 
hand,  is  unrealistic  for  non-uniform  microstructures,  e.g.  in  the  presence  of  clustering  of  heterogeneities  or 
microscopic  damage.  Even  with  a  uniform  phase  distribution  in  the  microstructure,  the  evolution  of  localized 
stresses,  strains  or  damage  path  can  violate  the  periodicity  conditions.  Problems  like  this  have  been  effectively 
tackled  by  multi-scale  modeling  methods  e.g.  in  [72,  29,  44,  67,  66,  81,  74,  73,  94,  108,  92].  Multi-scale 
analyses  methods  can  be  broadly  classified  into  two  classes.  The  first  is  known  as  ’’hierarchical  models” 
[29,  44,  94,  92]  in  which  information  is  passed  from  lower  to  higher  scales,  usually  in  the  form  of  material 
properties.  The  hierarchical  homogenization  models  assume  periodic  representative  volume  elements  (RVE) 
in  the  microstructure  and  uniformity  of  macroscopic  field  variables.  The  second  class,  known  as  “concurrent 
methods”  [67,  66,  82,  74,  73,  108],  implement  sub-structuring  and  simultaneously  solve  different  models  at 
regions  with  different  resolutions  or  scales. 

The  two-way  coupling  of  scales  enabled  in  the  concurrent  methods  is  suitable  for  problems  involving 
localization,  damage  and  failure.  Macroscopic  analysis,  using  bottom-up  homogenization  in  regions  of  rel¬ 
atively  benign  deformation,  enhances  the  efficiency  of  the  computational  analysis.  As  a  matter  of  fact,  it 
would  be  impossible  to  analyze  large  structural  regions  without  the  advantage  of  a  continuum  model  based 
macroscopic  analysis.  On  the  other  hand,  the  top-down  localization  process  cascading  down  to  the  mi¬ 
crostructure  in  critical  regions  of  localized  damage  or  instability  for  pure  microscopic  analysis,  is  necessary 
for  accurately  predicting  the  damage  path.  These  microscopic  computations,  depicting  the  real  microstruc¬ 
ture  are  often  complex  and  computationally  prohibitive.  Hence,  a  concurrent  setting  makes  such  analyses 
feasible,  provided  the  ”zoom-in”  regions  are  kept  to  a  minimum.  The  adaptive  multi-level  models,  promoted 
in  [67,  66,  82,  74,  73,  108],  are  attempts  to  achieve  this  objective,  with  the  adaptivity  motivated  from  phys¬ 
ical  and  mathematical  perspectives.  However,  there  is  a  paucity  of  such  studies  in  the  literature  involving 
material  nonlinearity  and  evolving  microstructural  damage.  In  their  previous  studies,  Ghosh  and  coworkers 
have  proposed  adaptive  multi-level  analysis  using  the  microstructural  Voronoi  cell  EEM  model  for  modeling 
elastic-plastic  composites  with  particle  cracking  and  porosities  in  [81],  and  for  elastic  composites  with  free 
edges  and  stress  singularities  in  [74,  73]. 

In  a  preceding  paper  [75],  the  authors  have  derived  and  computationally  modeled  an  anisotropic  con- 
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tinuum  damage  mechanics  (CDM)  model  for  unidirectional  fiber-reinforced  composites  undergoing  interfa¬ 
cial  debonding  from  by  using  homogenization  theory.  The  CDM  model  homogenizes  the  damage  incurred 
through  initiation  and  growth  of  interfacial  debonding  in  a  microstructural  RVE  with  nonuniform  distribu¬ 
tion  of  fibers.  Additionally,  arbitrary  loading  conditions  are  also  effectively  handled  by  this  model.  The 
present  paper  uses  this  CDM  model  of  [75]  in  an  adaptive  concurrent  multi-level  computational  model  to 
analyze  multi-scale  evolution  of  damage.  Damage  by  fiber-matrix  interface  debonding,  is  explicitly  modeled 
over  extended  microstructural  regions  at  critical  locations  [35,  53].  The  adaptive  model  addresses  issues  of 
efficiency  and  accuracy  through  considerations  of  physically-based  modeling  errors. 

The  adaptive  multi-level  model  consists  of  three  levels  of  hierarchy  viz.  level-0,  level- 1  and  level-2),  which 
evolve  in  sequence.  The  continuum  damage  model  developed  in  [75]  is  used  for  level-0  computations.  The 
level-1  domain  is  used  as  a  ‘swing  region’  to  establish  criteria  for  switching  from  macroscopic  to  microscopic 
calculations.  Physical  criteria  involving  variables  at  the  macroscopic  and  microstructural  RVE  levels,  trigger 
switching  from  pure  macroscopic  to  pure  microscopic  calculations,  i.e.  the  level  —  0  — >•  level  —  1  — >•  level  —  2. 
A  transition  layer  is  placed  between  the  level  —  1  and  microscopic  level  —  2  domains  for  smooth  transition 
from  one  scale  to  the  next.  All  computations  in  the  composite  microstructure  with  explicit  representations 
of  the  fiber  and  matrix  phases  are  done  with  the  Voronoi  cell  finite  element  model  or  VCEEM  [35,  53]. 
In  VCEEM,  debonding  at  the  fiber-matrix  interface  is  achieved  by  a  layer  of  cohesive  springs  [68].  Two 
numerical  examples  are  solved  in  this  paper  to  examine  the  effectiveness  of  the  multi-level  computational 
model  in  multi-scale  damage  analysis.  The  first  example  considers  a  small  region  of  a  fiber  matrix  composite 
microstructure  for  comparison  with  an  explicit  micromechanics  model.  The  second  set  of  problems  models 
a  double  lap  bonded  composite  joint  for  demonstrating  its  capability  in  handling  large  structural  problems. 

6.2  Levels  in  the  Multi-scale  Computational  Model 

The  multi-phase  composite  computational  domain  ^het  is  adaptively  decomposed  into  a  set  of  non-intersecting 
open  subdomains,  belonging  to  levels-0,  -1  and  -2  with  different  algorithmic  treatments,  i.e.  ^het  = 
;o  U  U  ^^2  U  ^tr  ■  The  different  levels  of  computational  hierarchy,  in  the  order  of  sequence  of  emer- 
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gence,  are  depicted  in  figure  (6.1)  and  discussed  briefly  here. 


LEVEL  2 


(a)  (b) 

Figure  6.1:  Schematic  of  the  two-way  coupled  concurrent  multi-level  model:  (a)  a  representative  volume 
element  (RVE)  for  a  non-uniformly  distributed  composite  microstructure  generated  by  tessellating  the  lo¬ 
cal  microstructure,  (b)  the  top-down  multi-level  model  showing  components  of  concurrent  coupling,  viz. 
continuum  level-0,  level-1  of  asymptotic  homogenization  and  level-2  of  micromechanical  analysis. 


6.2.1  Computational  Subdomain  Level-0  (Q/o  ) 

This  level  corresponds  to  regions  where  continuum  constitutive  laws  can  be  used  in  macroscopic  analysis. 
Macroscopic  field  variables  like  stresses  and  strains  in  Vlio  are  relatively  uniform  and  there  is  no  strong  non¬ 
periodicity  in  the  microstructure.  Hence,  microscopic  ‘statistical’  periodicity  in  the  RVE  is  assumed  to  be 
valid  in  this  level.  Scale  effects  are  negligible  and  it  is  possible  to  derive  effective  constitutive  relations  by 
volume  averaging  the  RVE  response  with  imposed  periodicity  conditions,  in  the  limit  that  the  RVE  tends  to 
zero  volume.  This  is  generally  the  starting  level  in  the  multi-scale  analysis  model,  as  long  as  RVE’s  can  be 
identified  for  the  computational  domain.  Macroscopic  analysis  with  the  continuum  constitutive  models  in 
level-0,  reduce  the  computing  effort  by  several  orders  of  magnitude  in  comparison  with  models  that  require 
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complete  microscopic  analysis  . 

For  undamaged  microstructures  with  linear  elastic  or  elastic-plastic  phases,  homogenized  anisotropic 
constitutive  laws  have  been  developed  by  the  authors  in  [82,  34].  In  the  case  of  microstructures  with  randomly 
evolving  microcracks  causing  diffused  damage,  the  homogenized  material  behavior  is  best  represented  by  a 
continuum  damage  mechanics  (CDM)  law.  An  anisotropic  CDM  model  with  a  fourth  order  damage  tensor 
has  been  developed  from  rigorous  micromechanical  analyses  in  [75].  The  general  form  of  CDM  models 
[49]  introduce  a  fictitious  effective  stress  acting  on  an  effective  resisting  area  (A),  which  is  caused  by 
reduction  of  the  original  resisting  area  A  due  to  material  degradation  from  the  presence  of  microcracks  and 
stress  concentration  in  the  vicinity  of  cracks.  In  [75],  the  effective  stress  is  related  to  the  actual  Cauchy 
stress  Sy  through  the  fourth  order  damage  effect  tensor  Mijki  as 

Sy  =  Mijki(D)'Ski  (6.1) 

where  Mijki  is  a  function  of  the  fourth  order  damage  tensor  D(=  Dijkiei  ®  ej  ®  ek  ®  ei).  The  hypothesis  of 
equivalent  elastic  energy  is  used  to  evaluate  Mijki  and  hence  establish  a  relation  between  the  damaged  and 
undamaged  stiffnesses  [18,  15,  106].  Equivalence  is  established  by  equating  the  elastic  energy  in  the  damaged 
state  to  that  in  a  hypothetical  undamaged  state  as 

iF(s,D)  =  isy(£:y«(D))-is«  =  lF(f:,0)  =  isy  (6.2) 

where  S  =  Syei®ej,  is  the  elastic  stiffness  tensor  in  the  undamaged  state  and  Eijki(D)  is  the  stiffness 
in  a  damaged  state.  From  equations  6.1  and  6.2,  the  relation  between  the  damaged  and  undamaged  stiffnesses 
is  established  as 

Eijkl  =  ^pqrsi^rskl)  (6-3) 

With  an  appropriate  assumption  of  a  function  for  Mijki,  equation  (6.3)  can  be  used  to  formulate  a  damage 
evolution  model  using  micromechanics  and  homogenization.  In  [75],  a  damage  evolution  surface  is  introduced 
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to  delineate  the  interface  between  damaged  and  undamaged  domains  in  the  strain  e-space  as 


F  = 


(6.4) 


Here  Wd  corresponds  to  the  dissipation  of  the  strain  energy  density  due  to  stiffness  degradation  for  constant 
strain  without  an  external  work  supply.  Also  called  the  degrading  dissipation  energy  (see  [43]),  it  is  an 
internal  variable  denoting  the  current  state  of  damage,  and  is  expressed  as: 


Wd 


f  1 


(6.5) 


Pijki  is  a  symmetric  negative-definite  fourth  order  tensor  that  will  be  expressed  as  a  function  of  the  strain 
tensor  aj,  a  is  a  scaling  parameter  and  k  is  a  function  of  Wd-  Assuming  associativity  rule  in  the  stiffness 
space,  the  evolution  of  the  fourth  order  secant  stiffness  is  obtained  as 


Pijki  —  ^ 


dF 


d{\eijeki) 


=  \Pi 


ijkl 


(6.6) 


Pijki  (e)  corresponds  to  the  direction  of  the  rate  of  stiffness  degradation  tensor  Eijki  -  For  a  composite  material 
with  interfacial  debonding,  the  direction  of  rate  of  stiffness  degradation  varies  with  increasing  damage  and 
hence  Pijki {e)  does  not  remain  a  constant  throughout  the  loading  process.  The  model  requires  the  evaluation 
of  K,  a  and  Ptjki  in  equation  (6.4).  These  are  determined  from  the  results  of  micromechanical  simulations  of 
a  RVE  with  periodic  boundary  conditions.  The  function  K{Wd)  is  evaluated  for  a  reference  loading  path  and 
all  other  strain  paths  are  scaled  with  respect  to  this  reference.  Upon  determination  of  the  maximum  value 
Wd  for  a  reference  loading  condition,  the  value  of  a  for  any  strain  path  can  be  obtained  by  simple  scaling. 
To  account  for  the  variation  of  Pijki  {e),  any  macroscopic  strain  evolution  path  is  discretized  into  a  finite  set 
of  points.  The  values  of  Pijki  are  explicitly  evaluated  at  these  points  from  RVE  based  simulations.  Values 
of  Pijki  for  any  arbitrary  macroscopic  strain  value  can  then  be  determined  by  interpolating  between  nodal 
values  using  shape  functions  of  a  3D  linear  hexahedral  element.  The  details  of  the  parameter  evaluation 
process  in  the  macroscopic  CDM  model  are  discussed  in  [75]. 
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6.2.2  Computational  Subdomain  Level-1 

Level-1  is  an  intermediate  computational  subdomain,  introduced  as  a  swing  region  for  establishing  criteria 
for  switching  from  macroscopic  level-0  regions  to  level-2  regions  of  pure  microscopic  computations.  The 
switching  criteria  are  based  on  analyses  of  the  macroscopic  problem,  as  well  as  of  the  microstructural  RVE 
problem.  The  asymptotic  homogenization  theory  is  used  for  this  level  to  decouple  the  set  of  governing 
equations  into  &  set  oi  (i)  homogenized  equations  representing  the  macroscopic  problem  corresponding  to  a 
length  scale  x,  and  (ii)  microscopic  equations  for  the  RVE  V(x),  represented  by  a  length  scale  y.  Details  of 
the  decoupled  macro-  and  micro-equations  are  given  in  the  appendix  section  6.7.1. 

Gradients  of  important  field  variables  are  evaluated  from  macroscopic  analysis  to  assess  the  deviation 
of  macroscopic  uniformity.  Such  gradients  may  be  the  effect  of  strong  microscopic  non-homogeneity  in  the 
form  of  highly  localized  stresses  and  strains  or  damage.  The  RVE-based  microscopic  analysis,  on  the  other 
hand,  provides  effective  criteria  to  estimate  departure  from  periodicity  conditions,  especially  in  the  event  of 
evolving  microstructural  damage.  The  adaptation  criteria  for  level  transitions  are  discussed  in  section  6.4. 
Two  sets  of  finite  element  problems  are  solved  for  the  level-1  subdomain  in  sequence,  viz., 

1.  Macroscopic  analysis:  Incremental  macroscopic  analysis  of  the  computational  domain  is  performed  using 
the  CDM  model  to  evaluate  macroscopic  variables  e.g.  stresses  and  strains  due  to  the  increments  in 
applied  loads. 

2.  Microstructural  RVE  analysis:  This  is  a  post-processing  operation  in  which  microstructural  analysis 
of  the  RVE  is  conducted  for  each  integration  point  of  the  macroscopic  elements.  The  strain  field  e^, 
obtained  from  macroscopic  analysis  with  the  CDM  model,  is  imposed  on  the  RVE  as  an  external  driver, 
together  with  periodic  boundary  conditions  on  the  boundary  of  the  RVE  as  shown  in  figure  (6.1)a. 
Microscopic  stresses  ,  strains  £ij  and  other  variables  are  computed  in  this  post-processing  stage  for 
each  RVE. 

Remark  1:  The  macroscopic  computations  of  level-0  and  level-1  elements  are  performed  with  the  conventional 
displacement-based  finite  element  method,  while  all  microscopic  calculations  in  the  RVE  of  level-1  elements 
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are  performed  using  the  Voronoi  cell  FEM  [80,  35,  53]. 

Remark  2:  Computational  models  in  the  macroscopic  level-0  and  level-1  subdomains  are  refined  adaptively  by 
selective  h-  or  h-p  strategies.  ‘Error’  and  convergence  criteria  for  this  refinement  have  been  discussed  in  [74]. 
Local  enrichment  through  successive  mesh  refinement  or  enhancement,  serves  a  dual  purpose  in  the  multi¬ 
level  computational  strategy.  The  first  goal  is  to  identify  regions  of  high  discretization  ‘error’  and  improve 
convergence  through  mesh  enhancement.  The  second  is  to  identify  regions  of  high  modeling  error  and  zoom 
in  on  these  regions  to  create  higher  resolution.  These  regions  are  generally  characterized  by  large  gradients 
and  localization  of  macroscopic  variables.  Element  refinement  in  these  regions  is  helpful  for  reducing  the 
length-scale  difference  between  macroscopic  elements  in  the  homogenized  domain  and  microscopic  regions 
with  explicit  representation  of  heterogeneities. 

6.2.3  Computational  Subdomain  Level-2  (Q/2  ) 

The  level-2  subdomain  of  pure  microscopic  analysis  emerges  from  level-1  elements  in  regions  characterized 
by  (a)  departure  from  macroscopic  uniformity,  e.g.  regions  of  localization  or  fracture,  and  (b)  significant 
microstructural  non-uniformities  manifested  by  e.g.  growth  of  localized  damage.  Prior  to  transition  to  level- 
2  elements,  a  high  spatial  resolution  is  reached  in  the  macroscopic  mesh,  resulting  in  small  elements,  by 
h-  or  hp-  refinement.  The  successive  refinement  process  stops  when  a  certain  element  size  is  achieved  and 
subsequently  the  model  changes  from  macroscopic  to  pure  microscopic.  A  scale  ratio  SR  is  chosen  a-priori 
to  ascertain  this  element  size.  Depending  on  the  choice  of  SR  =  oY  focal  j  microscopic 

model  in  any  given  level-2  element  can  encompass  large  portions  of  the  microstructure  with  many  discrete 
heterogeneities.  The  elements  are  constructed  by  filling  with  the  exact  microstructure  at  that  location, 

as  outlined  in  the  following  steps  and  shown  in  figure  6.2. 

•  Use  appropriate  adaptation  criteria  to  determine  if  a  level-1  element  needs  to  switch  to  level-2  element. 

•  Identify  a  region  in  the  microstructure  ^Imicro  that  is  located  in  the  same  region  as  the  level- 2  element. 
^micro  should  extend  beyond  the  element  boundary  by  approximately  two  fiber  lengths. 

•  Tessellate  the  local  microstructure  to  generate  a  mesh  of  Voronoi  cell  elements  as  shown  in  figure  (6.3). 
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•  Carve  out  the  microstructural  region  of  the  level- 2  element  from  the  local  microstructure  ^micro-  This 
procedure  will  result  in  dissecting  some  of  the  fibers  on  the  boundary.  When  this  happens,  additional 
nodes  are  generated  on  the  Voronoi  cell  boundary  at  locations  where  the  fiber  surface  and  Voronoi 
cell  edges  intersect  the  boundary  of  the  level- 2  element.  The  dissected  conjugate  pieces  of  a  fiber 
belonging  to  two  contiguous  level-2  elements  are  joined  together  when  the  two  contiguous  elements 
share  a  common  edge. 


o 


□  Transition  element 
O—Q  Special  interface  layer 


□  Level-0/ 1  Element 

□  Level-2  element 
•  Level-0/1  nodes 

O  Level-0/1  nodes  at  the  transition  interface 
■  VCFEM  nodes  on  level-2/transition  boundary 
X  VCEEM  internal  nodes 
Y  Transition  element  nodes  at  the  interface 


(a)  (b) 

Figure  6.2:  (a)  Process  of  carving  out  level-2  element  microstructure  (b)  Interface  constraints  between  level- 
0/level-l  and  tr  elements 

Requirement  of  high-resolution  micromechanical  models  in  these  elements  entails  prohibitively  large  com¬ 
putations  using  conventional  finite  element  methods.  The  microstructure-based  Voronoi  cell  FEM  [35,  53,  80] 
is  particularly  effective  for  modeling  level-2  elements  because  of  its  efficiency  in  modeling  large  heterogeneous 
regions  [35,  53,  80,  81,  74].  Each  Voronoi  cell  with  embedded  heterogeneities  (particle,  fiber,  void,  crack  etc.) 
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represents  the  region  of  contiguity  for  the  heterogeneity,  and  is  treated  as  an  element  in  VCFEM.  VCFEM 
elements  can  be  considerably  larger  than  conventional  EEM  elements  and  incorporate  a  special  hybrid  EEM 
formulation.  Incorporation  of  known  functional  forms  from  analytical  micromechanics  substantially  enhances 
its  convergence.  A  schematic  diagram  of  Voronoi  cell  elements  is  shown  in  figure  (6.3).  A  high  level  of  ac¬ 
curacy  has  been  achieved  with  VCEEM  for  modeling  problems  with  microstructural  damage  by  particle 
cracking  [36]  and  fiber- matrix  interfacial  debonding  [35,  53].  Eor  debonding  simulation,  imperfect  interfaces 
are  represented  by  the  cohesive  zone  model  [68].  Displacement  degrees  of  freedom  on  the  fiber- matrix  in¬ 
terface  are  constrained  by  the  cohesive  zone  models  as  discussed  in  section  6.3.  VCEEM  has  been  shown 
to  be  significantly  more  efficient  than  commercial  displacement  based  EE  packages  for  modeling  complex 
microstructures  with  evolving  damage. 


Figure  6.3:  A  typical  level-2  element  containing  an  aggregate  of  microstructural  Voronoi  cell  elements  with 
relevant  notations. 
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6.2.4  Scale  Transition  Subdomain  (Qtr) 


The  interface  between  the  level-0  or  level-1  elements  and  the  level-2  elements  with  explicit  representation 
of  the  heterogeneous  microscopic  domain,  needs  a  special  treatment  to  facilitate  smooth  transition  of  scales 
across  the  element  boundaries.  A  layer  of  transition  elements  {Etr  G  ^Itr)  is  sandwiched  between  these 
elements,  where  (fljr)  is  the  transition  subdomain  as  shown  in  figure  (6.2)b.  The  E^r  elements  are  essentially 
level-2  elements  with  compatibility  and  traction  continuity  constraints  imposed  at  the  interface  with  level- 
0/level-l  elements.  It  is  assumed  that  layers  of  Etr  elements  are  located  beyond  the  critical  hot-spots, 
at  which  homogenization  fails.  Hence,  the  homogenized  laws  are  sufficient  at  their  interfaces  with  level- 
1/level-O  elements.  A  weak  form  of  the  interface  displacement  continuity  is  incorporated  through  the  use  of 
Lagrange  multipliers  on  this  interface  [74,  73].  This  results  in  a  weak  satisfaction  of  the  interface  displacement 
compatibility  and  avoids  the  occurrence  of  spurious  forces  that  arise  if  the  displacements  are  strongly  coupled. 


6.3  Coupling  Different  Levels  in  the  Concurrent  Multi-Scale  Al¬ 
gorithm 

The  concurrent  multi-scale  analysis  requires  that  all  levels  be  coupled  for  simultaneously  solving  for  variables 
in  the  different  computational  subdomains.  Consequently,  the  global  stiffness  matrix  and  load  vectors  are 
derived  for  the  entire  computational  domain  {^Ihet  =  U  Ha  U  ^li2  U  The  corresponding  domain 

boundary  is  delineated  as  Thet  =  {Lo  U  Ta  U  ri2}  where  T^o  =  dUio  n  Thet]  La  =  dUn  n  Thet]  ^12  = 
d^li2  n  Tftet-  Let  Tint  =  dTlii  n  dTltr  delineate  the  boundary  between  the  level- 1  and  transition  elements, 
where  the  displacement  continuity  is  satisfied  using  Lagrange  multipliers.  The  incremental  form  of  the 
equation  of  principle  of  virtual  work  equation  for  Tlhet  at  the  end  of  an  increment,  can  be  written  as  the  sum 
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of  contributions  from  each  individual  domain,  as 


I  (Si,-  +  ASi,-)  — ^  dn-  {U  +  Mi)5uf  dr 

Ow  JTio 


f  f 

+  /  (Si,-  +  ASi,-)— ^  df]  -  {ti  +  Ati)6u‘l  dr 
Jnn  Jrn 


r 

(ctj,-  -I-  Actj,-)  — — ^  dU  - 
n./  9xj  j 

f  {ti  +  Ati)dwf  dr 

r*. 

1  {Oij  +  Act*,-)  '  dn  1 

^  {ti  -\-  Ati)5uf  dr 
ri2 

f  +  AAf  )(^;i  +  A^;i  -  wf  -  Awf  )dr  =  0 


)dr 


(6.7) 


The  prefix  A  symbolizes  increments  of  the  respective  variables  in  the  incremental  solution  process.  The 
superscripts  10,  11,  12,  tr  correspond  to  association  -with  the  respective  level,  while  the  (/)  sign  refers  to 
variables  that  could  belong  to  either  level.  Si,-  are  the  components  homogenized  macroscopic  stresses  ob¬ 
tained  from  the  CDM  constitutive  model  for  rtio  and  rta .  The  applied  tractions  ti  are  at  traction  boundaries 
of  the  respective  domains.  The  boxed  parts  in  equation  (6.7)  correspond  to  contributions  from  level-2  and 
transition  computational  subdomains  that  are  generated  from  VCFEM  solutions  of  the  microstructural  re¬ 
gions.  Displacement  components  rtf,  and  uf  are  on  the  boundaries  of  elements  coinciding  with  the 
boundaries  of  the  riio,  rin,  and  0^2  subdomains.  An  intermediate  segment  Tjnj  is  added  at  the  interface 
between  the  level-1  and  tr  elements,  as  shown  in  figure  6.2.  On  these  segments,  displacement  components  Vi 
are  interpolated  with  any  order  polynomial  functions,  independent  of  the  interpolations  on  or 

Even  for  highly  nonhomogeneous  displacements,  high  order  interpolations  on  the  intermediate  segment  are 
able  to  smoothen  the  transition  between  levels.  This  has  been  demonstrated  for  problems  without  damage 
through  numerical  examples  in  [74].  The  last  two  terms  in  equation  (6.7)  use  Lagrange  multipliers  to  fa¬ 
cilitate  incorporation  of  a  weak  form  of  the  interfacial  displacement  continuity  on  Ejnj.  and  are 

vector  columns  of  Lagrange  multipliers  belonging  to  domains  Dio/a  and  ritr  respectively  at  rmt-  The  Euler’s 
equations,  obtained  from  setting  the  coefficients  of  Svi,  and  to  zero  respectively  in  the  principle 
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of  virtual  work  (6.7),  are 


^m/a  ^  ^  ^  =  _^xt  +  AAf)  =  -{a^j  +  Aa.jY^nj 

(ui  +  =  (ui  +  Auj)*’’  =  {vi  +  Avi)  (6.8) 

where  n*  is  the  unit  normal  vector  and  and  XY  correspond  to  the  interfacial  traction  components  on 

dUiQ/ii  and  dUtr  respectively.  The  displacements  Vi  and  the  Lagrange  multipliers  xf^‘^  and  xf^*^  on  the 
intermediate  boundary  segment  are  interpolated  from  nodal  values  using  suitably  assumed  shape  functions 
as: 

{v}  =  [Lint]{(lint}  ,  =  [I/;^iO/il]{A;o/a}  !  (6.9) 

The  displacements  uf  and  in  each  level-0  and  level-1  elements  are  interpolated  by  the  standard  or 
hierarchical  Legendre  polynomials  based  shape  functions  as: 

{u}™  =  [N,o]{q.o}  =  [Nfo  NO]  < 

{u}'i  =  [Na]{qa}  =  K  NO]  < 

As  shown  in  figure  6.2,  the  generalized  displacements  in  the  level-0  and  level-1  elements  are  subdivided  into 
two  classes:  (i)  those  at  nodal  points,  which  interface  with  transition  elements,  and  (ii)  those  at  all  other 
nodes.  Generally,  only  level-1  elements  will  interface  with  transition  elements  because  of  the  sequence  of 
introduction  of  computational  levels.  The  generalized  displacements  qfo/a  corresponds  to  the  nodal  degrees 
of  freedom  in  level-O/level-1  elements  at  the  interface  with  transition  elements,  while  q^/;i  correspond  to 
the  remaining  degrees  of  freedom  in  these  elements.  The  solution  of  the  algebraic  form  of  equation  (6.7)  is 
obtained  using  the  Newton  Raphson  iterative  solver.  Setting  up  the  tangent  stiffness  matrix  requires  consis¬ 
tent  linearization  by  taking  directional  derivative  of  equation  (6.7)  along  incremental  displacement  vectors 


q/o 


q/i 


(6.10) 
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Au  and  Av,  and  the  Lagrange  multipliers  AA.  For  the  i— th  iteration  in  the  solution  of  the  incremental 


variables,  assembled  matrix  equations  derived  from  equation  (6.7)  has  the  following  structure. 
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As  explained  before,  the  superscript  I  represents  quantities  associated  with  nodal  points  at  the  interface  with 
transition  elements  while  superscript  O  indicate  association  with  nodes  at  other  regions.  The  two  notations 
in  the  superscript  separated  by  comma,  represents  the  node  coupling  effect.  For  example,  the  superscript 
7,0  corresponds  to  the  coupling  between  the  non-interface  and  interface  nodes.  The  stiffness  submatrices 
[K;o/a]  and  sub-vector  {F;o/a}  correspond  to  those  for  the  level-0  and  level-1  elements  and  are  expressed  as 
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(6.12) 


The  subscripts  {m,n)  correspond  to  the  degrees  of  freedom  while  {a,  (3)  correspond  to  the  node  numbers  in 
the  element.  These  matrices  and  vectors  are  further  divided  based  on  the  classification  of  the  I  and  O  nodes. 
The  coupling  between  the  level- 0 /level- 1  and  tr  elements  is  achieved  through  the  [P]  and  [Q]  matrices,  which 


134 


may  be  expressed  as 


JTir^t 

^Ptr)moin^ 

JTi„t 

iQ  10  /  ll')  manP  —  /  (lL)  ma  (Ij^l0/ll)nl3dT 
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(6.13) 


Contributions  to  the  load  vector  {F }  due  to  coupling  between  level- 0 /level- 1  and  tr  elements  are  given  as 


JTi 

(iL)„(A'°/''  +  AX‘°/‘^)mdT  -  [  (Lf„J„(A'2/*’-  +  AA'2/*’-)„dr 

nt  Jri„t 

{^^XlO/ll)ma  =  — 

(6.14) 

{^^Xl2/tr)ma  —  ~ 

hir^t 

Finally,  the  stiffness  [Ki2/tr]  and  the  load  vector  {Fi2/tr}  for  level-2  and  tr  elements  are  obtained  by  VCFEM 
calculations  followed  by  static  condensation  to  represent  the  virtual  work  in  terms  of  the  boundary  terms 
only. 


6.3.1  Modified  Voronoi  Cell  FEM  Formulation  for  a  RVE  in  Level-1  Elements 

Details  of  the  Voronoi  Cell  FEM  are  provided  in  [80,  35,  53]  and  have  been  summarized  in  the  appendix  section 
6.7.2.  As  discussed  in  section  6.2.2,  the  post-processing  phase  for  level-1  elements  require  the  evaluation  of 
different  variables  in  the  RVE  from  known  values  of  macroscopic  strains.  A  small  variant  of  the  formulation  in 
equation  (6.40)  enables  this  execution.  The  energy  functional  for  a  RVE  (Y)  with  V-periodic  displacements 
and  V-anti-periodic  tractions  on  the  boundary,  and  imposed  macroscopic  strain  (e^j  -I- Ae^),  may  be  written 
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as 
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(6.16) 


The  boxed  term  corresponds  to  the  additional  energy  due  to  the  imposed  macroscopic  strain  field  on  the 
RVE  region  Y.  The  Euler-Lagrange  equations  corresponding  to  this  energy  functional  are: 


£ij(x,y)  +  A£y(x,y)  =  Sijki{(Tij  +  Act^)  =  (ey(x)  +  Aeij(x)) 
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Ui  is  E-periodic  and  cr^rij  is  E-anti-periodic  on  dY^ 


(6.17) 

(6.18) 


The  corresponding  weak  form  of  the  element  kinematic  relation  is  written  in  a  matrix  equation  form  as 
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or  in  a  condensed  form 


[H'=]{;a  +  A^a}  =  +  Aq}  -  {Rf }  (6.20) 

This  relation  is  then  substituted  in  equation  (6.44)  for  obtaining  the  RVE  based  solutions.  It  should  be 
noted  that  displacement  periodicity  is  imposed  on  the  RVE  boundary  for  solving  this  problem. 

6.4  Criteria  for  Adaptive  Mesh  Refinement  and  Level  Transitions 

In  the  application  of  the  multi-level  model,  the  following  criteria  are  used  for  mesh-refinement  and  level 
transitions  due  to  discretization  and  modeling  error  respectively.  Many  of  these  adaptation  criteria  are 
based  on  the  physics  of  the  problem  in  consideration,  since  rigorous  mathematical  error  bounds  are  scarce 
(or  even  non-existent)  for  these  nonlinear  problems.  Consequently  they  are  nonunique  and  other  indicators 
may  be  used  if  appropriate. 

6.4.1  Refinement  of  Level-0  and  Level-1  Meshes  by  h-Adaptation 

The  computational  models  in  level-0  and  level-1  subdomains  are  enriched  by  h-adaptation  based  mesh  re¬ 
finement  to  reduce  discretization  ‘error’.  The  /i— adaptation  procedure  subdivides  candidate  macroscopic 
elements  into  smaller  elements  to  reduce  a  suitably  chosen  error.  It  is  necessary  to  impose  boundary  dis¬ 
placement  compatibility  constraint  conditions  between  contiguous  divided  and  undivided  elements  in  this 
method  [82].  This  local  mesh  enrichment  is  intended  to  reduce  discretization  error  and  to  identify  regions 
of  modeling  error  by  zooming  in  on  localization  regions  with  evolving  gradients.  Eor  CDM  based  evolving 
problems,  an  adaptation  criterion  is  formulated  in  this  paper  in  terms  of  the  jump  in  traction  across  adjacent 


137 


element  boundaries  that  signifies  local  stress  gradients.  The  condition  is  stated  as: 


Refine  element  ‘e’,  if  the  traction  jump  error  across  the  element 
satisfies  the  condition:  >  Ci  *  , 


where 


f  ptj\2 

E%  =  and  {E^^)l 
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lan^  {m?  +  [[Ty]?)ddn 

fdOe 


(6.21) 


Here  NE  is  the  total  number  of  level-0  and  level- 1  elements  in  the  entire  computational  domain,  T^^Ty  are 
the  components  of  element  boundary  tractions  in  the  x  and  y  directions  and  [[.]]  is  the  jump  operator  across 
element  boundary  dVle-  A  factor  Ci  (<  1)  has  been  chosen  from  numerical  experiments. 


6.4.2  Criteria  for  Switching  from  Level-0  to  Level-1  Elements 

Level-0  to  level-1  element  transition  takes  place  according  to  criteria  signaling  departure  from  conditions  of 
the  homogenizability  that  are  based  on  macroscopic  variables  in  the  continuum  model  of  level-0  elements. 
The  degrading  dissipation  energy  Wa  in  the  CDM  model  is  a  strong  indicator  of  localized  damage  evolution. 
Consequently,  a  criterion  is  formulated  as: 


Switch  element  ’e’  from  level-0  to  level-1  if  : 

*  iWa)e  >  C2  *  E^J:,  *  {Wd)ma.  (6.22) 


where  ® 


is  the  norm  of  the  local  gradient  of  the  degrading  dissipation  energy  {Waje,  expressed  as: 
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Emax  i®  111®  maximum  value  of  for  all  elements  and  C2  (<  1)  is  a  prescribed  factor.  The  criterion  (6.22) 
is  helpful  for  seeking  out  regions  with  high  gradients  of  Wa  in  regions  of  high  Wa  itself.  In  a  previous  paper 
by  the  authors  [74],  the  gradient  was  expressed  in  terms  of  the  maximum  difference  in  the  damage  for 
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neighboring  elements  as  =  Max|(h7d)e  —  (hrd)adjacentl-  A  more  accurate  definition  of  the  local  gradient 
is  adopted  in  the  present  work,  using  the  Zienkiewicz-Zhu  (ZZ)  gradient  patch  recovery  method  [107].  In 
this  method,  interpolation  of  Wd  is  assumed  in  the  form  of  a  polynomial  over  a  patch  of  elements  adjoining 
a  nodal  point  in  a  level-0  element.  The  least  square  minimization  process  leads  to  the  local  matrix  equation 


e(xi ,  X2)]^[Ne(xi ,  X2)]{a}  =  Yl[^e{xuX2)f{Wd)e{xi,X2)  (6.24) 

e=l  e=l 

where  [Ne(a:i,X2)]  is  a  matrix  containing  polynomial  interpolation  terms  and  ne  is  the  number  of  elements 
in  the  patch.  The  equation  (6.24)  is  solved  for  the  coefficients  {a}.  The  gradients  of  Wd  in  each  element  are 
calculated  from  the  nodal  values  using  element  shape  functions  as 


dWd 

dxi 


dWd 

dX2 


ON, 

dx2 


iWd)c 


(6.25) 


6.4.3  Criteria  for  Switching  from  Level- 1  to  Level-2  Elements 

For  elements  in  which  macroscopic  nonuniformity  has  been  established  according  to  equation  (6.22),  de¬ 
parture  from  RVE  periodicity  is  taken  as  an  indicator  for  activating  a  switch  from  level-1  to  level-2.  The 
switching  criterion  is  developed  in  terms  of  evolving  variables,  e.g.  the  averaged  strain  at  the  fiber-matrix 
interface  in  the  local  microstructural  RVE.  The  averaged  strain  is  stated  as: 


"  7 — r  ^  (6-26) 

Judiia  JudQc  JudOa  ludQc 

where  the  integral  is  evaluated  over  all  the  fiber-matrix  interfaces  in  the  RVE.  The  jump  in  displacement 
across  the  fiber- matrix  interface  with  a  normal  rij  is  denoted  by  \uj\.  Eor  perfect  interfaces  \uj\  will  be  zero. 
Thus,  Dij  corresponds  to  the  contributions  to  macroscopic  strain  due  to  damage  only,  and  Dij  =  0  in  the 
absence  of  damage.  Departure  from  periodicity  will  result  in  a  significantly  different  averaged  strain  Dij 
in  response  to  different  conditions  on  the  boundary  of  the  microstructural  region.  Eor  example,  let 
correspond  to  the  solution  of  a  boundary  value  problem  of  the  local  microstructure  included  in  a  level-2 
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element  (see  figure  6.2)  subject  to  boundary  displacements  that  have  been  obtained  from  the  macroscopic 
level-0/1  analysis.  The  scale  of  the  microstructure  is  relevant  in  this  analysis  since  periodicity  is  not  imposed 
on  the  boundary.  On  the  other  hand,  let  be  from  the  solution  of  a  boundary  value  problem  of  the 

local  RVE  with  imposed  macroscopic  strains  and  subjected  to  periodic  boundary  displacements  constraints. 
The  difference  in  these  two  strains  for  a  level-1  element  e  may  be  quantified  as 


=  max(|Ti"’'^ 
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(6.27) 


For  evaluating  in  a  given  step  of  the  incremental  solution,  only  the  increments  in  the  present  step 
are  calculated  by  the  level-1  macroscopic  displacement  boundary  conditions.  It  is  assumed  that  the  RVE- 
based  solution  is  valid  all  the  way  upto  (but  excluding)  the  present  step.  The  departure  from  periodicity  is 
measured  in  terms  of  the  difference  in  averaged  strains  and  hence  the  criterion. 


Switch  element  ’e’  from  level-1  to  level-2  if  : 

(6.28) 

where  is  the  maximum  value  of  for  all  the  level- 1  elements  in  the  computational  domain. 

Remark  Once  the  regions  of  level-2  and  transition  elements  have  been  identified,  it  is  important  to  update 
the  local  micromechanical  states  of  stress,  strain  and  damage  to  the  current  state.  This  step  should  precede 
the  coupled  concurrent  analysis.  For  this  analysis,  the  history  of  the  macroscopic  displacement  solution  on 
the  level-0 /level- 1  element  boundary  prior  to  the  switch  is  utilized.  The  local  micromechanical  (VCFEM) 
boundary  value  problem  for  the  level-2  element  is  incrementally  solved  from  the  beginning  to  obtain  the 
history  of  stresses,  strains  and  damage  in  the  microstructure  from  the  macroscopic  boundary  displacement 
history. 
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6.5  Numerical  Examples  with  the  Adaptive  Multilevel  Model 


Two  sets  of  numerical  examples  are  solved  to  study  the  effectiveness  of  the  multi-level  computational  model 
for  composite  materials. 

6.5.1  Multi-level  Model  vs.  Micromechanical  Analysis 

This  example  is  aimed  at  understanding  the  effectiveness  of  the  multi-level  model  in  analyzing  a  nonuniform 
composite  microstructure  by  comparing  its  predictions  with  those  by  pure  micromechanical  analysis.  It  is 
computationally  intensive  to  conduct  pure  micromechanical  analysis  with  evolving  damage  for  very  large  mi- 
crostructural  regions.  Consequently  a  computational  domain  with  a  small  population  of  fibers,  as  shown  in 
the  optical  micrograph  of  figure  6.4(a),  is  considered.  The  micrograph  is  for  a  polymer  matrix  composite  with 
a  random  dispersion  of  uniaxial  fibers.  The  dimensions  of  the  micrograph  analyzed  are  100/<m  x  70.09/<m, 
containing  264  circular  fibers  of  diameter  1.645/<m  each,  corresponding  to  a  volume  fraction  of  32%.  Though 
the  domain  may  not  be  adequate  for  a  clear  separation  between  continuum  and  micromechanical  regions 
(since  relatively  large  regions  are  needed  to  materialize  the  RVE) ,  the  results  of  this  example  are  enough  to 
show  the  effectiveness  of  the  overall  framework. 

The  optical  micrograph  is  mapped  onto  a  simulated  microstructure  with  circular  fibers  that  is  tessellated 
into  a  mesh  of  264  Voronoi  cell  elements,  as  shown  in  figure  6.4(b).  The  constituent  materials  in  the  composite 
system  are  an  epoxy  resin  matrix,  stainless  steel  reinforcing  fibers  and  a  very  thin  film  of  freekote  (<  0.1/<m) 
at  the  fiber-matrix  interface.  The  freekote  imparts  a  weak  strength  to  the  steel-epoxy  interface,  which  allows 
a  stable  growth  of  the  debond  crack  for  experimental  observation.  The  experimental  methods  of  material  and 
interface  characterization  have  been  discussed  in  [35].  Both  the  matrix  and  fiber  materials  are  characterized 
by  isotropic  elasticity  properties.  The  matrix  material  has  a  Young’s  modulus,  Egpoxy  =  4.6  GPa  and 
Poisson’s  ratio,  Vepoxy  =  0.4,  while  the  fiber  material  has  a  Young’s  modulus,  Egteei  =  210  GPa  and 
Poisson’s  ratio,  Vsteei  =  0.3.  A  bilinear  cohesive  law  described  in  [53,  68]  is  used  in  this  analysis  for  modeling 
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(a)  (b) 

Figure  6.4:  (a)  Optical  micrograph  of  a  steel  fiber-epoxy  matrix  composite  with  264  fibers  (b)  the  simulated 
computational  model  with  a  Voronoi  cell  mesh 


the  fiber-matrix  interface.  In  this  model,  the  normal  and  tangential  tractions  are  given  as 
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where  t  is  a  bilinear  function  of  the  interfacial  separation  as 


(6.29) 
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The  unloading  behavior  in  the  hardening  region  is  linear  following  the  loading  path.  In  the  softening  region, 
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the  unloading  proceeds  along  a  different  linear  path  from  the  current  position  to  the  origin  with  a  reduced 
stiffness,  for  which  the  t  —  6  relation  is 


t  = 


^max  ^max  '■ 
^max 


6c  <6 

max  <  6e  and  6  <  6„ 


(6.31) 


It  is  expected  that  the  degrading  dissipation  energy  Wd  in  the  macroscopic  CDM  model  depends  on  the 
cohesive  parameters  in  the  microstructural  debonding  model.  A  square  RVE  with  a  single  circular  fiber  is 
simulated  for  interfacial  debonding  with  three  different  sets  of  cohesive  parameters,  as  shown  in  the  inset 
of  figure  6.5.  The  cohesive  energies  are  the  same  for  all  cases.  However  in  one  case,  the  critical  separation 
length  6e  is  increased  while  in  the  other,  the  corresponding  peak  stress  a^ax  is  reduced.  The  figure  6.5  infers 
that  while  6e  has  a  small  influence  on  Wd,  the  effect  of  amax  is  certainly  significant,  at  least  in  the  early 
stages  of  straining. 


The  cohesive  parameters  used  in  this  paper  are:  amax  =  0.005  GPa,  dc  =  5.1  x  10^®  m  and  6e  = 
3.1  X  10^^  m.  The  microstructure  is  loaded  in  tension  in  the  horizontal  direction  with  a  displacement  of 
0.1/<m  that  is  uniformly  increased  in  20  equal  increments,  corresponding  to  a  total  strain  of  en  =  0.1%.  The 
displacement  boundary  condition  is  imposed  along  the  right  edge,  as  shown  in  figure  6.4(b). 

Micromechanical  analysis  by  VCFEM 

The  pure  micromechanical  VCEEM  solution  using  the  mesh  of  figure  6.4(b)  has  been  presented  in  [53]  and 
are  used  here  as  reference  solutions  for  the  multi-scale  simulation.  Figure  6.6(a)  shows  the  contour  plot  of 
microscopic  stress  a^x  at  the  final  step  of  the  micromechanical  simulation  with  a  depiction  of  interfacial 
debonding.  The  right  side  of  the  microstructure  shows  significant  concentrated  damage  with  this  load.  The 
debonding  initiates  at  the  top  and  percolates  to  the  bottom  of  the  microstructure  along  a  narrow  band. 

Multi-Scale  Analysis  with  the  Multi-Level  Model 

Multi-scale  analysis  is  performed  by  the  concurrent  multi-level  computational  model  and  the  results  are 
compared  with  those  from  the  micromechanical  VCFEM  analysis.  For  the  multi-level  model,  the  entire  com- 
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2.5e-05 


Figure  6.5:  The  degrading  dissipation  energy  Wd  as  a  function  of  strain,  evaluated  for  different  cohesive  zone 
parameters  in  the  bilinear  cohesive  law 

putational  region  of  264  fibers  is  first  divided  into  9  macroscopic  finite  elements  as  shown  in  figure  6.7(a).  For 
evaluating  the  homogenized  constitutive  properties  for  each  of  element,  statistically  equivalent  representative 
volume  element  or  SERVE  for  the  microstructure  underlying  each  macroscopic  element  is  first  identified. 
Various  statistical  methods  have  been  used  to  determine  the  size  scale  of  the  RVE  and  the  number  of  inclu¬ 
sions  contained  in  it  [70,  37,  86,  98].  Rigorous  methods  of  evaluating  statistically  equivalent  representative 
volume  elements  by  a  combination  of  statistical  methods  and  micromechanical  analyses  have  been  conducted 
by  the  first  author  in  [83,  90].  However,  since  the  number  of  fibers  in  the  micrograph  is  limited  in  this  ex¬ 
ercise,  a  simpler  assumption  is  made.  The  SERVE  for  each  element  is  assumed  to  consist  of  all  the  fibers 
belonging  to  that  element.  Eor  example,  to  generate  the  SERVE  for  an  element  window  in  the  micrograph 
of  figure  6.4(b),  all  fibers  whose  centers  are  located  within  this  window  are  first  identified  as  constituents 
of  the  RVE.  This  is  shown  by  the  aggregate  of  black  fibers  in  figure  6.8(a).  The  homogenization  method, 
discussed  in  sections  6.2.2  and  6.3,  requires  a  periodic  distribution  of  the  RVE  and  this  is  achieved  by  locally 
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repeating  the  arrangement  of  fibers  in  both  the  xi  and  X2  directions  for  a  period  length  in  figure  6.8(a).  This 
means  that  for  each  fiber  identified  in  the  element,  at  (xi ,  X2),  four  identical  fibers  are  placed  at  the  locations 
(xi  ± Xi,  X2),  (xi,  X2  ±  X2)  where  (Xi,  X2)  are  periods  in  the  two  directions.  The  period  lengths  Xi,  X2  are 
selected  such  that  the  volume  fraction  of  RVE  matches  that  of  the  local  microstructure.  Finally,  the  domain 
is  tessellated  into  a  network  of  Voronoi  cells  as  shown  in  figure  6.8(b)  Tessellation  provides  a  natural  way 
of  creating  periodic  SERVE  boundary.  For  non-uniform  fiber  arrangements,  the  SERVE  boundary  consists 
of  non-straight  line  edges.  The  nodes  on  this  SERVE  boundary  are  periodic,  i.e.  for  every  boundary  node 
a  periodic  pair  can  be  identified  on  the  boundary  at  a  distance  of  one  period  along  each  of  the  coordinate 
directions.  In  figure  6.8(b),  the  node  pairs  are  identified  as  AA,  BB  etc.  The  number  of  fibers  and  their 
distribution  in  the  SERVE  of  each  macroscopic  element  is  shown  in  figure  6.7(a). 

Since  the  number  of  elements  in  this  exercise  is  very  small  (only  9),  level-0  simulations  with  the  CDM 
model  is  bypassed  in  the  multi-level  analysis.  All  elements  are  level-1  at  the  start  of  the  multi-level  simulation. 
Switch  to  level-2  elements  is  made  in  accordance  with  equations  (6.27)  and  (6.28)  with  C2  =  0.2.  However 
the  terms  for  each  element  in  equation  (6.27)  are  replaced  by  the  difference  in  RVE  based 

averaged  strains  between  adjacent  elements  Also,  as  opposed  to  an  entire  macroscopic 

element,  a  single  layer  of  transition  Voronoi  cell  elements  is  included  between  the  level-1  and  level-2  elements. 
In  figure  6.7(b)  the  Voronoi  elements  containing  the  grey  fibers  constitute  the  transition  layer,  while  those 
containing  the  black  fibers  belong  to  level-2.  An  interface  segment  Tint  is  inserted  between  the  transition  and 
level-1  elements  at  a  distance  Ltr/12  from  the  right  edge.  Convergence  properties  of  the  multi-level  model 
are  studied  by  considering  two  cases  with  =  o.35  and  =  0.45.  This  is  achieved  by  changing  the 

size  of  the  initial  level-1  elements. 

As  depicted  in  figure  6.7(b),  only  three  elements  (3,  6  and  9)  at  the  right  side  of  the  initial  mesh  switch 
from  level-1  to  level-2.  A  comparison  of  results  by  (a)  VCFEM  based  micromechanical  analyses  (all  level-2 
elements)  ,  (b)  homogenization  based  macroscopic  analysis  (all  level- 1  elements),  and  (c)  concurrent  multi¬ 
level  analysis  {level-1  and  level-2  elements)  for  =  o.35  and  0.45  is  made.  Contour  plots  of  cth  (GPa) 
showing  interfacial  debonding  at  the  end  of  the  simulation  are  shown  for  the  concurrent  multi-scale  analysis 
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in  figure  6.6(b,c).  The  discrepancy  in  the  damage  path  predicted  by  the  micromechanical  analysis  and  the 
multi-level  analysis  reduces  sharply  with  increasing  value.  This  can  be  attributed  to  the  fact,  that 

the  damage  path  is  very  sensitive  to  the  macro-micro  interface  conditions.  Since  the  sample  size  is  small 
and  there  is  no  real  periodicity  in  the  microstructure,  the  proximity  of  the  level-1  boundary  to  the  damage 
localization  zone  alters  the  local  boundary  conditions.  However  as  this  distance  is  increased,  the  microscopic 
stress  distribution,  debonding  pattern  and  damage  zone  replicates  the  real  event  observed  in  micromechanical 
analysis.  The  distribution  of  the  micromechanical  stresses  cth  ,  generated  by  pure  micromechanical  and  multi¬ 
level  analyses,  are  plotted  along  a  line  through  the  middle  of  micrograph  in  figure  6.9.  The  micromechanical 
stresses  show  only  minor  oscillations  about  an  averaged  value  of  the  0.005  GPa  in  the  region  to  the  left 
of  the  level- 1 -level- 2  interface.  In  the  region  to  the  right,  where  damage  is  predominant,  there  is  clearly  a 
convergence  of  the  stresses  with  increasing  value. 

The  macroscopic  or  averaged  stress-strain  response  for  element  1  (always  level-1)  and  element  9  (changes 
levels)  are  plotted  in  figures  6.10.  For  the  micromechanical  problems  with  debonding,  the  volume  averaged 
stresses  and  strains  are  evaluated  by  averaging  the  local  fields  over  the  microscopic  domain  as: 

^  /  (^ij{xi,X2)dQ,  and  ^ij  =  ^  [  ^ij{xi,X2)d^  -  Dij  (6.32) 

where  Dij  is  the  strain  jump  defined  in  equation  (6.26).  The  results  for  all  the  models  are  in  good  agreement 
for  the  element  1,  where  there  is  no  significant  microstructural  damage.  The  small  difference  is  due  to  the 
periodicity  constraints  imposed  on  the  microstructure.  Also  there  is  a  difference  between  the  results  of  case 
1:  =  0.35  and  case  2:  =  0.45,  due  to  the  interface  conditions  at  Vint-  However,  as  is  expected 

the  results  are  quite  different  for  element  9,  where  significant  damage  is  observed  in  figure  6.6.  The  level-1 
analysis  shows  significant  deviation  from  the  micromechanical  analysis  due  to  imposed  periodicity  in  the 
damage  zone.  Once  again,  the  results  improve  significantly  with  increasing  ratio. 
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6.5.2  A  Composite  Double  Lap  Joint  with  Microstructural  Debonding 

Adhesive  bonded  joints  are  considered  as  preferred  alternatives  to  fasteners  for  joining  structural  components 
due  to  their  light  weight.  However,  adhesively  bonded  structures  consisting  of  different  materials,  can  induce 
high  stresses  near  the  interface  leading  to  failure  initiation  by  interfacial  debonding.  A  double-lap  bonded 
joint  with  boron-epoxy  composites  as  adherents,  is  analyzed  in  this  example.  An  adhesive  shown  as  ABCD 
in  figure  6.11(a)  is  used  to  bond  the  two  composite  materials.  Only  a  quarter  of  the  joint  is  modeled  from 
considerations  of  symmetry  in  boundary  and  loading  conditions.  For  boundary  conditions,  the  displacement 
component  ui  is  set  to  zero  along  the  face  X2  =  0  implying  symmetry  about  the  Xi  axis.  The  displacement 
components  ui  and  M2  along  the  face  xi  =  8/i  are  set  to  zero  corresponding  to  a  fixed  edge.  A  tensile 
displacement  ui  is  applied  on  the  face  of  the  lower  ply  at  xi  =  0.  Both  plies  above  and  below  the  adhesive 
are  made  of  unidirectional  boron  fiber-  epoxy  matrix  composite  materials.  The  fibers  are  uniformly  arranged 
in  a  square  array  in  the  microstructure,  implying  a  square  unit  cell  with  a  single  circular  fiber.  The  epoxy 
matrix  has  a  Youngs  modulus  E  =  4.6  GPa  and  Poisson’s  ratio  v  =  0.4,  while  boron  fibers  have  a  Youngs 
modulus  E  =  210  GPa  and  Poisson’s  ratio  v  =  0.3.  The  material  properties  of  the  isotropic  adhesive  are: 
Young’s  modulus  E  =  3.45  GPa  and  Poisson’s  ratio  v  =  0.34.  The  bilinear  cohesive  law  parameters  for  the 
matrix-fiber  interface  are:  CTtoox  =  0.02  GPa,  6c  =  5.0  X  10  ®  m  and  Sg  =  20.0  x  10  ^  m. 

Multi-level  analysis  for  model  with  450  fibers 

In  this  model,  the  top  ply  above  the  adhesive  consists  of  10  rows  of  fiber,  while  the  bottom  row  consists  of 
5  rows  resulting  in  a  total  of  450  fibers.  The  microstructural  volume  fraction  of  fibers  is  Vf  =  20%.  The 
applied  displacement  on  the  face  at  Xi  =  0,  is  uniformly  increased  from  zero  to  ui  =  1.2  x  10^®/i  in  15 
uniform  increments.  The  number  of  fibers  is  kept  low  in  this  example,  such  that  a  micromechanical  analysis 
can  be  easily  done  for  this  example  with  a  mesh  of  450  Voronoi  elements,  each  of  which  is  a  square  unit  cell. 
The  micromechanics  solutions  are  used  as  a  reference  to  determine  the  accuracy  of  multi-scale  simulations. 
Three  different  approaches  are  used  to  solve  this  problem.  They  are:  (a)  a  macroscopic  model  using  the 
continuum  damage  model  for  constitutive  behavior,  (b)  a  detailed  micromechanical  VCFEM  analysis,  and 
(c)  a  multi-level  model  for  multi-scale  analysis.  The  starting  mesh  in  the  multi-level  model  of  the  bonded 
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joint  consists  of  a  uniform  grid  of  470  QUAD4  elements  for  macroscopic  analysis  as  shown  in  figure  6.11. 
The  constitutive  relation  for  each  element  is  a  fourth  order  anisotropic  CDM  model  that  has  been  devel¬ 
oped  for  this  unit  cell  with  interfacial  cohesive  zone  in  [75].  Figure  6.12(a)  shows  the  contour  of  degrading 
dissipation  energy  Wd,  at  the  final  stage  of  loading  by  a  pure  CDM  based  macroscopic  analysis.  Damage 
initiates  near  the  bottom  left  corner  A  of  the  adhesive  joint  and  propagates  downwards  to  span  the  entire 
region  on  the  left  of  point  A.  Level-0  — >•  level-1  transition  in  the  multi-level  analysis  is  performed  using 
equation  (6.22)  and  level-1  — >•  level-2  transition  uses  equation  (6.28)  with  factors  C2  =  0.5  and  C3  =  0.1. 
The  gradient  of  the  energy  ^ (^^)^  +  (^7)^  the  final  loading  stage,  used  in  equation  (6.22),  is  shown 
in  figure  6.12(b).  The  corresponding  evolution  of  various  levels  in  the  multi-scale  model  is  depicted  in  fig¬ 
ure  6.13  at  two  different  loading  stages.  There  are  7  level-1  elements  at  87%  of  the  final  loading.  At  the  final 
load  increment,  the  multi-level  mesh  consists  of  446  level-0  elements,  0  level-1  elements,  14  level-2  elements 
and  10  transition  elements.  All  level- 2  elements  emerge  in  critical  the  regions  where  both  the  gradient  and 
intensity  of  Wd  are  high  in  the  macroscopic  analysis.  Figure  6.14(a,b)  depict  the  contours  of  microscopic 
stress  CTii  and  the  regions  of  debonding  obtained  by  pure  micromechanical  and  the  multi-level  models.  The 
results  of  the  multi-level  model  are  in  excellent  agreement  with  the  micromechanical  analysis,  both  with 
respect  to  debonding  regions  and  evolving  variables.  The  maximum  error  in  cth  is  around  1%.  The  excellent 
agreement  is  further  corroborated  in  the  plot  of  cth  along  the  vertical  line  through  the  microstructure  in 
figure  6.9.  Figure  6.10(a,b)  plot  the  macroscopic  (averaged)  Sn  —  en  curve  obtained  from  (a)  macroscopic 
CDM-based  analysis,  (b)  micromechanical  analysis  and  (c)  multi-scale  analysis  with  the  multi-level  model 
at  two  different  locations,  PI  and  P2  shown  in  figure  6.11(b).  At  P2,  where  the  damage  and  its  gradient 
are  low,  solutions  by  the  CDM  model  and  micromechanics  are  in  relatively  good  agreement.  At  this  point, 
the  multi-scale  model  uses  the  CDM  constitutive  law.  However,  the  CDM  results  are  quite  different  from 
the  other  two  at  PI,  a  hotspot  where  the  damage  and  its  gradient  are  high.  It  is  assuring  to  note  that  the 
multi-level  model  matches  the  micromechanics  results  quite  well  at  this  point. 

The  computational  efficiency  of  the  multi-level  model  is  examined  by  a  comparison  of  the  CPU  time 
on  a  IA32  computer  cluster  for  the  different  models.  The  computations  are  carried  out  in  a  serial  manner 
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using  a  single  processor.  The  results  are  tabulated  in  table  6.1.  Although  the  macroscopic  CDM  analysis  is 
faster,  it  can  lead  to  significant  errors.  The  complete  level-1  solution  is  even  slower  than  the  micromechanics 
solution,  since  it  solves  the  RVE  problem  in  every  element.  Accurate  analysis  with  the  multi-level  model  is  at 
least  7  times  faster  than  the  complete  micromechanics  and  level-1  solutions  for  this  problem.  The  efficiency 
increases  rapidly  with  increasing  number  of  fibers  in  the  analysis. 


Model 

Level-0 

Level- 1 

Micromechanics  (Level-2) 

Multi-scale 

Time  in  seconds 

71 

300330 

300310 

42260 

Table  6.1:  CPU  time  on  a  IA32  cluster  to  solve  the  double  lap  joint  model  by  various  methods. 

Multi-level  analysis  for  model  with  192000  fibers 

This  is  a  more  realistic  model  of  the  composite  joint  with  a  large  number  of  fibers,  to  realize  the  potential 
of  the  multi-level  model.  The  top  ply  consists  of  160  rows  of  fiber,  while  the  bottom  row  consists  of  80  rows 
resulting  in  a  total  of  192000  fibers.  The  geometric  and  material  parameters  are  the  same  as  in  the  previous 
example,  except  for  the  special  cases  mentioned.  A  pure  micromechanical  analysis  is  not  conducted  due 
to  the  large  number  of  fibers.  The  problem  is  analyzed  by  (a)  a  macroscopic  model  by  CDM  and  (b)  the 
multi-level  model.  The  multi-level  analysis  activates  all  three  types  of  adaptation: 

•  Refinement  of  level-0  elements  by  h-adaptation  in  accordance  with  equation  (6.21),  for  Ci  =  0.7. 

•  Transition  from  level-0  to  level-1  elements  in  accordance  with  equation  (6.22),  with  C2  =  0.5. 

•  Transition  from  level- 1  to  level-2  elements  in  accordance  with  equation  (6.28),  with  C3  =0.1. 

The  effects  of  variation  of  cohesive  zone  parameters  and  the  effect  of  volume  fraction  are  studied.  The  unit 
cells  considered  in  this  example  have  2  volume  fractions:  (i)  Vf  =  20%.  and  (ii)  Vf  =  40%.  Three  different 
cases  with  different  parameters  in  the  bilinear  cohesive  law  are  considered. 

•  Cl:  Same  cohesive  parameters  as  in  section  6.5.2. 

•  C2:  Umax  and  5e  are  the  same  as  in  section  6.5.2.  However,  5c  is  4  times  that  in  case  Cl.  This  reflects 
the  same  cohesive  energy  with  a  smaller  ascending  slope. 
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•  C3:  Umax  is  reduced  by  half  and  6^  is  doubled.  Hence  the  cohesive  energy  is  the  same  as  Cl  with  a 
smaller  peak  stress.  Also  5c  is  the  same  as  that  in  Cl. 

The  starting  mesh  has  470  level-0  elements.  For  Vf  =  40%  and  case  Cl,  the  final  mesh  has  1688  level-0 
elements,  24  level-2  elements  and  33  transition  elements  as  shown  in  figure  6.15(a).  Figure  6.15(b)  illustrates 
the  corresponding  microscopic  stress  distribution  and  debonding  in  the  level-2  regions  near  the  hotspot  at 
A.  The  macroscopic  (averaged)  stress-strain  plots  are  shown  for  two  points  in  the  composite  joint:  (a)  near 
the  critical  point  A  and  (b)  at  a  non  critical  point  B  are  shown  in  figure  6.16.  The  predictions  of  the  CDM 
model  agree  with  the  multi-level  model  at  the  point  B.  However,  the  stress  predictions  by  the  CDM  model 
are  considerably  higher  than  those  by  the  multi-level  model  at  A,  where  damage  is  very  localized  and  the 
periodicity  condition  imposed  by  the  CDM  model  is  unrealistic. 

The  effect  of  Vf  on  the  damage  evolution  near  the  corner  PI  is  seen  in  figure  6.17  for  the  case  Cl. 
A  significantly  higher  Wa  is  observed  for  the  higher  volume  fraction,  which  increases  with  evolving  strain. 
Figure  6.18  shows  the  distribution  of  Wa  at  the  end  of  the  analysis  for  the  different  cohesive  parameters. 
Intense  damage  localization  takes  place  near  the  junction  A  in  the  bond  (see  figure  6.11(b).  Damage  starts 
from  this  location  and  propagates  down  and  left  towards  the  edge  of  the  applied  loading.  Damage  localization 
is  the  strongest  for  the  case  Cl,  and  propagates  almost  vertically  down  in  a  narrow  zone.  It  is  in  these  regions, 
that  scale  transition  to  level- 2  occurs.  The  damage  distribution  in  the  remaining  parts  of  the  composite  joint 
is  rather  low  and  uniform.  Moving  the  peak  stress  in  case  C2  with  a  lower  traction-displacement  slope  results 
in  a  more  diffused  damage  region  and  the  damage  seem  to  spread  more  in  the  region  to  the  left  of  point 
A.  The  damage  localization  reduces  for  the  case  C3  with  lower  peak  stress  and  the  damage  is  more  evenly 
distributed.  For  Vf  =  20%,  the  damaged  regions  are  less  localized. 

6.6  Conclusions 

An  adaptive  concurrent  multi-level  computational  model  is  developed  in  this  paper  for  multi-scale  analysis 
and  prediction  of  damage  in  fiber  reinforced  composite  materials.  Microstructural  damage  is  manifested  by 
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fiber-matrix  interfacial  debonding  in  this  paper.  Microstructural  damage  mechanisms  leading  to  complete 
failure  are  more  complex  than  the  singular  mode  of  damage  considered  in  this  paper.  The  authors  are 
currently  working  towards  this  goal,  where  interfacial  debonds  bifurcate  into  the  matrix  and  eventually 
coalesce  to  cause  a  continuous  fracture  path.  A  step  forward  in  this  direction  can  be  seen  in  a  recent  paper 
on  the  growth  and  coalescence  of  multiple  cohesive  cracks  [54].  However,  the  intent  of  the  present  paper  is  to 
create  a  framework  for  the  multi-scale  coupling  so  that  more  complex  damage  mechanisms  may  eventually 
be  incorporated.  Hence  interfacial  debonding  is  deemed  sufficient  for  this  purpose. 

The  multi-level  model  invokes  two-way  coupling  of  scales,  viz.  a  bottom-up  coupling  with  homogenization 
at  lower  scales  to  introduce  reduced  order  continuum  models  and  a  top-down  coupling  at  critical  hotspots 
to  transcend  scales  for  following  the  microstructural  damage  evolution.  The  bottom-up  coupling  results  in  a 
continuum  damage  mechanics  (CDM)  model  developed  in  a  preceding  paper  [75].  Three  levels  of  hierarchy, 
with  different  resolutions,  evolve  in  this  model  with  adaptation.  Adaptive  capabilities  enable  effective  domain 
decomposition  in  the  evolving  problem  with  damage,  keeping  a  balance  between  computational  efficiency  and 
accuracy.  Macroscopic  analysis  is  done  with  the  CDM  model  of  [75]  for  high  efficiency.  Pure  micromechanical 
analysis  is  computationally  exhaustive  and  the  adaptive  methodology  optimally  reduces  this  region  to  a 
minimum.  The  Voronoi  cell  finite  element  model  [35,  53]  is  effectively  utilized  for  efficient  micromechanical 
analysis  of  extended  microstructural  regions.  The  numerical  examples  establish  the  accuracy  and  efficiency 
aspects  of  the  model,  as  well  as  demonstrate  its  capability  in  handling  problems  involving  damage  in  large 
composite  domains.  Overall  this  work  lays  an  effective  foundation  for  solving  multi-scale  problems  involving 
localization,  damage  and  crack  evolution  that  may  be  impossible  to  achieve  using  any  single  scale  model. 

6.7  APPENDIX 

6.7.1  Microscopic  and  Macroscopic  Equations  in  Computational  Subdomain 
Level- 1 

Any  function  /  in  the  RVE  is  assumed  to  be  T-periodic,  i.e.  /(x,  y)  =  /(x,  y  -I-  kY)  Vfc  =  1,  2,  •  •  • .  Peri¬ 
odicity  conditions  are  used  on  the  RVE  boundary  to  decouple  the  set  of  equations  at  different  levels  as: 
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Microscopic  equations 


eij(x,y)  =  ey(x)  +  -( 


l^dui{y)  ,  duj{y) 


2^  dyj 


+ 


dyi 


o-y(x,y) 

9CTij(x,y) 

dyj 


eu  (x)  +  ^ki  (*)  (Kinematics) 

/y\ 

EijkiSmkSni  +  e^^j.(x)^^^e„„(x)  (Constitutive) 


=  0  (Equilibrium) 


(6.33) 


Macroscopic  equations 


'(x)  “  ]^  /  +  ^kL~^)  dYemn  =  Efj^^emni^)  (Constitutive) 


^Sij(x) 

dxi 


+  /*  =  0  (Equilibrium) 


(6.34) 


In  the  above  equations  Ui  is  a  E-periodic  displacement  function  and  (Tjj(x,  y)  is  the  stress  field  in  the  RVE 
respectively,  while  Sjj(x)  and  e^-  are  the  homogenized  stress  and  strain  tensors.  Eij^i  and  Efjf,i  correspond 
to  microscopic  and  homogenized  anisotropic  elasticity  tensor  respectively.  The  details  of  the  derivation  of 
equations  (6.33)  and  (6.34)are  discussed  in  [81,  74]. 


6.7.2  The  Voronoi  Cell  Finite  Element  Model 

A  typical  level-2/  transition  element  consisting  of  microstructural  Voronoi  cell  elements  is  shown  in  figure 
6.3.  The  following  assumptions  are  made  in  the  formulation  of  each  VCEEM  element. 

•  Stress  fields  a/-  in  the  matrix  phase  and  in  the  inclusion  phase  fig  of  each  Voronoi  cell  element 
fie  are  independent  and  equilibrated.  The  stress  interpolations  in  each  phase  are  expressed  as  : 


{(7  ™}  =  [P™]{;a™}  in  fi„  and  {(T  “}  =  [P=]{;a=}  in  fie 


(6.35) 


where  the  matrices  [P™A]  are  obtained  from  assumed  stress  functions  like  the  Airy’s  stress  function 
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and  {P  ™/c|  aj-g  unknown  coefficients  to  be  solved. 


•  Compatible  displacement  field  are  assumed  on  each  Voronoi  cell  element  boundary  dQ.e  and  inter¬ 
polated  as: 

{u-}  =  [L"]{q"}  (6.36) 

•  Compatible  displacement  fields  u™  and  u?  are  assumed  on  the  matrix  and  inclusion  sides  of  the  matrix- 
inclusion  interface  respectively,  and  are  interpolated  as  : 

{u™}  =  [L=]{q™}  on  SO™  and  {u=}  =  [L=]{q=}  on  (6.37) 

In  an  incremental  formulation,  the  potential  energy  functional  for  each  element  is  expressed  in  terms  of  the 
incremented  stresses  and  displacements  as: 


ne(at^,  Au?)  =  -  [  ,  Aa^)dU 

-  [  AB%atj,Aatj)dn  +  [  (a^  +  Aa^)(e^  +  Ae^)di} 

JQc 

r  r  r'^'n  ~'^n~ 

+  /  (af,.  +  Aaf,.)(e?,-  +  Aey)dn  -  /  /  T™d(u™  -  Odd^ 

JQc  J  dQc 

ru^  -\-Au^  —Uf  —  Auf 


-/  / 

J  dO-c  uV 


^d{u^  -  ui)ddn  -  [  {ti  +  Ati){i 

JTtrr. 


?  +  A<)dr 


(6.38) 


Here  B  =  \Sijki(Jij(Jki  is  the  complementary  energy  density  and  AB  =  ^SijkiAuijAuki  +  SijkiAuijUki-  The 
strain  fields  and  are  in  the  matrix  and  inclusion  phases  respectively  of  each  Voronoi  element,  t  is  the 
prescribed  traction  on  the  boundary  Ttm-  The  prefix  A  corresponds  to  increments  and  subscripts  n  and  t 
correspond  to  the  normal  and  tangential  directions  at  the  matrix-inclusion  interface.  The  two  terms  on  the 
matrix-inclusion  interface  DO.™ IdQ.^  provide  the  work  done  by  the  interfacial  tractions  T™  =  T™n™  -|-T/"t™ 
due  to  interfacial  separation  (u™  —  u'^).  T™  and  T™  are  the  normal  and  tangential  components  that  are 
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described  by  cohesive  laws  in  [35,  53].  Using  divergence  theorem,  the  potential  energy  can  be  written  as: 


ne  =  -  /  dn-  [  dn 

~  f  dU  —  f  S^jf.iahAa'^j  dfi 

+  f  (a["  +  Aa^)n]iu^i  +  Au^)ddn  -  f  (a[”  +  Aa[”)n=«  +  AuT)ddn 

JdQe  Jan™ 


(6.39) 


+  (atj  +  Aatj)n%u<i  +  Au<i)ddn  -  /  T™d«  -  <)d50 

JdQ^.  JdUcJu^-u^ 


p  pU^  i-ZAXAj  ~'^t 

-  /  /  Trd{uT  -  uDdon 

J dQ.c 

-[  iU  +  Ati)iu^i  +  Au^i)dT 


(6.40) 


Here  n®  and  n®  are  the  outward  normal  on  9fle  and  9flc  respectively.  The  integration  over  the  incremental 
displacements  at  the  interface  9flc  is  conducted  by  the  backward  Euler  method.  The  total  potential  energy 
functional  for  each  level-2  or  tr  element  containing  N^c  Voronoi  cell  elements  as  shown  in  figure  6.3  is 

iV„c 

n'2/*’'  =  ^ne  (6.41) 

e=l 


Substituting  stress  interpolations  (6.35)  and  displacement  interpolations  (6.36,6.37)  in  equation  (6.40)  and 
setting  variations  with  respect  to  the  stress  coefficients  A/3™,  A/3®  respectively  to  zero  results  in  the  weak 
form  of  the  element  kinematic  relation. 


[P™]^[S™][P™]dO 

[0] 


[0] 

/o  [Pl^[Sl[Plci« 


^a®  +  A^a® 


>  = 


[P™]^[n®][L®]dafl  -  [P™]^[n®][L®]dafl 


[0] 


[0] 


[0] 

/9QjP®]^[n®][L®]dafl 


q®  +  Aq® 
q™  +  Aq” 
q®  +  Aq® 


or  in  a  condensed  form 


[H®]{;a  +  A;a}  =  [G®]{q  +  Aq} 


(6.42) 
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The  weak  forms  of  the  global  traction  continuity  conditions  are  subsequently  solved  by  setting  the  variation 
of  the  total  energy  function  in  equation  (6.41)  with  respect  to  Aq,  Aq™  and  Aq*^  to  zero.  This  results  in 
the  weak  form  of  the  traction  reciprocity  conditions  as: 


iV„e 

E 

e=l 


[0] 

[0] 

[0] 


/B’^  +  A/B 

jB^  +  AjB 


>  = 


J\vc 

E 


< 


e=l 


ijLr{i+Ai}d(} 

-  Jqq^  ({n‘=}T™(M„  +  Aun,ut  +  Aut)  +  {t=}T/"(M„  +  Am„,  ut  +  Aut))  ddO,  * 

-  ({n®}T™(M„  +  Au„,  Ut  +  Aut)  +  +  Am„,  ut  +  Aut))  ddO. 

(6.43) 


or  in  a  condensed  form 


Nvc 


Nvc 


Y^[Gr{i3+m=Ei^^} 


(6.44) 


Substituting  (6.42)  in  (6.44)  yields: 


f^[Gl^[Hl-MGl{q  + Aq}  = 


e=l 


e=l 


(6.45) 


In  an  iterative  solution  of  equation  (6.45),  its  linearized  form  for  the  i— th  iteration  is  given  as: 


iV„e  iV„e  iV„e 

^[Gl^[H'=]-i[Gl)HAq}*  =  -  £[Gl^[ff]-l[Gl{q  +  Aq}* 

e=l  e=l  e=l 

or  in  a  condensed  form 


(6.46) 


[K]*{Aq}*  =  {AF*""}* 


(6.47) 


In  order  to  incorporate  this  relation  in  the  linearized  form  of  the  principle  of  virtual  work  of  equation  (6.11), 
it  should  be  noted  that  the  displacement  vector  Ui2/tr  on  the  boundary  of  level- 2  and  transition  element  is 
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a  subset  of  all  the  VCFEM  displacement  fields,  i.e.  =  u®  U  u™  U  u®.  Consequently,  the  displacement 
field  u^®  can  be  separated  into  two  categories,  viz.  (i)  nodal  points  at  the  boundary  of  the  level-2 

or  transition  elements  shown  is  figure  6.3,  and  (ii)  u*"*  on  all  the  internal  nodes.  The  stiffness  matrix  and 
the  load  vector  of  the  ensemble  of  all  Voronoi  cell  elements  belonging  to  a  level-2  or  transition  element  can 
therefore  be  partitioned  as 


J^12/tr,12/tr  J^12/tr,int 

t 

Aq>2/tr 

t 

^F>2/tr 

< 

>  =  < 

j^int,l2/tr 

Aq‘"* 

<  j 

AF’"* 

<  > 

Static  condensation  of  the  internal  degrees  of  freedom  yields 


=  {AF“*}* 


(6.48) 


(6.49) 


These  stiffness  matrices  and  load  vectors  are  then  used  in  global  assembly  process  of  equation  (6.11). 
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Figure  6.6:  Contour  plot  of  cth  showing  interfacial  debonding  at  the  end  of  the  simulation,  for:  (a)  pure 
micromechanical  analysis,  (b)  analysis  by  multi-scale  model  with  a  smaller  level-2  region  =  0.35),  and 

(c)  analysis  by  multi-scale  model  with  a  larger  level-2  region  =  0.45). 


Figure  6.7:  Computational  mesh  for  the  computational  domain:  (a)  Macroscopic  mesh  with  different  RVE 
in  every  element,  (b)  Multi-level  model  with  the  interface  between  macroscopic  and  microscopic  VCFE 
elements. 


SERVE  boundary 

D  1  C 


Level-1  element  boundary 


Eigure  6.8:  (a)  A  periodic  microstructure  containing  the  tessellated  RVE  (fibers  in  black),  (b)  placement  of 
the  RVE  in  the  level-1  element  showing  periodic  nodes  on  the  boundary. 


Figure  6.9:  Comparison  of  microscopic  stress  an  by  different  methods,  plotted  along  a  line  through  the 
middle  of  microstructure 
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(a) 


(b) 

Figure  6.10:  Comparison  of  macroscopic  (volume  averaged)  Sn  —  en  curves  by  different  methods  of  analysis 
at  (a)  macroscopic  element  1,  and  (b)  macroscopic  element  9. 
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Figure  6.11:  (a)  Schematic  diagram  of  a  composite  double  lap  joint  showing  dimensions  and  boundary 
conditions,  (b)  the  level-0  computational  mesh. 
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(b) 

Figure  6.12:  Contour  plot  of  (a)  degrading  dissipation  energy  Wa  and  (b)  its  gradient 
at  the  final  loading  stage. 
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INCREMENT  13/15 


TOTAL  ELEMENT  470 

□  LEVEL  0  =  461 

□  LEVEL  1  =  7 

□  TRANSITION  ELEMENTS  =  5 
■  LEVEL  2  =  4 


(a) 

INCREMENT  15/15 
TOTAL  ELEMENT  470 

□  LEVEL  0  =  446 

□  LEVEL  1  =  0 

□  TRANSITION  ELEMENTS  =  10 
■  LEVEL  2=  14 


(b) 

Figure  6.13:  Evolution  of  the  multi-level  computational  model  with  level  transition  (a)  at  87  %  loading,  and 
(b)  at  the  final  loading  stage. 
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(b) 

Figure  6.14:  Level  2  microscopic  VCFEM  elements  near  the  corner  A  showing  microscopic  stress  distribution 
(GPa)  and  interfacial  debonding  at  the  end  of  the  analysis  by:  (a)  pure  micromechanical  analysis  (b)  multi¬ 
scale  analysis. 
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INCREMENT  18/18 
TOTAL  ELEMENT  1745 

□  LEVEL  0  =  1688 

■  LEVEL  1  =  0 

□  TRANSITION  ELEMENTS  =  33 

■  LEVEL  2=  24 


(a) 


Max. 


Min. 


(b) 


Figure  6.15:  (a)  Evolved  multi-level  model  and  mesh  at  the  final  load  step,  (b)  Microscopic  stress  distribution 
and  interfacial  debonding  at  the  end  of  analysis  for  location  near  corner  A,  for  Vf  =  40%  and  case  Cl. 
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0.05 


(a) 


(b) 

Figure  6.16:  Macroscopic  averaged  stress-strain  (Sn  —  en)  plot  at  two  locations  in  the  double  lap  joint:  (a) 
critical  region  A  and  (b)  non-critical  region  B. 
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Damage  Work 


Figure  6.17:  Degrading  dissipation  energy  evolution  near  the  corner  A  of  the  double  lap  joint  for  Vf  =  20% 
and  Vf  =  40%,  and  case  Cl. 
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Figure  6.18:  Distribution  of  Wd  with  Vf  =  40%  and  different  cohesive  parameters:  (a)  case  Cl,  (b)  case  C2, 
and  (c)  case  C3,  at  the  end  of  loading. 
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